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Summary 

Work  was  conducted  on  two  projects  related  to  control  of  shear  flows. 

In  the  first,  two-dimensional  simulations  were  performed  to  study  the  effects  of  resolution  on 
the  extraction  of  velocity  data  from  measurements  of  scalar  fields  in  incompressible  flows.  That 
work,  with  Bonnie  Carpenter,  has  been  published  in  Physics  of  Fluids,  8,  2447-2459  (1996)  and 
demonstrates  "proof-of  concept"  for  the  methodology  described  in  an  earlier  Physics  of  Fluids 
paper,  also  supported  by  AFOSR.  A  reprint  is  included  in  this  report. 

In  the  second,  David  Cotrell  and  I  have  developed  a  general-purpose  algorithm  for  computing 
all  solutions  of  an  arbitrary  system  of  multivariate  polynomial  equations.  That  work  constitutes  the 
basis  of  a  technique  for  determining  the  coefficients  in  low-order  models  of  scalar  transport  in 
turbulent  shear  flows.  The  polynomial  root-finder  work  (preprint  included  in  this  report)  will  be 
submitted  to  the  SIAM  Journal  on  Scientific  Computing.  The  application  to  low-order  model 
construction  is  described  in  a  paper  (with  Ari  Glezer,  Martin  Brooke,  and  their  students)  to  be 
submitted  to  Journal  of  Fluid  Mechanics). 

1.  Global  Scalar  Velocimetry 

This  work  constituted  the  bulk  of  Bonnie  Carpenter's  MS  Thesis.  A  summary  is  provided 
here. 

Under  earlier  AFOSR  support,  we  developed  a  technique  for  determining  n-dimensional 
(/i=2,3)  velocity  fields  from  measurement  of  n-l  passive  or  reactive  scalars  has  recently  been 
developed  (Pearlstein  and  Carpenter,  Phys.  Fluids,  7,  754-763,  1995).  That  method  utilizes  n 
linear,  first-order,  uncoupled  hyperbolic  equations  for  n  velocity  components,  derived  from 
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transport  equations  for  n-\  scalars  and  conservation  of  mass.  The  velocity  is  determined  by 
integration  along  the  characteristic  curves  defined  by  the  hyperbolic  equations.  In  a  second  paper 
{Phys.  Fluids,  8,  2447-2459,  1996),  supported  by  Grant  F49620-95-1-0297,  we  presented  the 
results  of  a  computational  proof-of-concept  study  for  a  steady,  two-dimensional,  diverging  channel 
flow.  We  considered  the  effects  of  grid  resolution  (and  the  shape  of  its  elements)  on  which  the 
scalar  is  known  and  its  derivatives  are  approximated,  the  integration  step  size  used  to  calculate  the 
velocity  components  along  characteristic  curves,  and  the  effects  of  multiplicative  noise  introduced 
into  the  scalar  field.  The  results  show  that  extraction  of  the  velocity  by  integration  of  the  equations 
along  characteristics  is  stable,  and  that  the  techniques  proposed  (ibid.)  for  removal  of  singularities 
are  effective.  For  steady  flows,  noise  in  the  scalar  field  measurements  can  be  dealt  with  by 
extracting  the  velocity  from  the  mean  of  several  noisy  scalar  fields. 

2.  Low-Order  Models  of  Scalar  Transport  in  Turbulent  Flows 

2.7.  Fixed-Point  Algorithm  for  Solution  of  Multivariate  Polynomial  Equation  System 
We  have  developed  a  globalization  of  the  Krawczyk  algorithm  to  compute  all  real  isolated 
solutions  of  systems  of  N  real  polynomial  equations.  This  is  accomplished  by  transforming  the 
original  system  to  an  augmented  system  in  7?^^* ,  in  which  each  of  the  first  TV  variables  lies  in  the 
interval  [-1,1],  and  the  (TV+l)-th  variable  lies  in  [0,~).  The  domain  of  this  latter  variable  can  be 
divided  into  two  intervals  [0,1]  and  [1,°°),  the  second  of  which  is  mapped  to  [0,1].  Thus,  the 
entire  domain  7?^'^'  can  be  examined  by  considering  the  finite  domain  [-1,1]^  X  [0,1]  for  each  of 
two  systems.  One  of  the  augmented  systems  can  have  one  or  more  solutions  not  corresponding  to 
finite  solutions  of  the  original  system.  We  use  techniques  from  algebraic  geometry  to  transform 
the  (7V+l)-th  dimensional  system  so  that  spurious  solutions  are  excluded,  thus  restricting  the 
solutions  to  those  corresponding  to  finite  solutions  of  the  original  system.  The  algorithm  requires 
bounds  for  multivariate  polynomials  over  a  finite  domain.  The  best  of  three  bounding  methods 
considered  uses  an  interval  extension  of  the  system,  which  is  stored  as  a  rooted,  ordered  tree,  and 
is  equivalent  to  an  TV-dimensional  Homer  scheme  that  takes  advantage  of  polynomial  sparsity. 
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We  have  extended  the  approach  to  deal  with  systems  having  nonsimple  solutions  (i.e.,  with 
multiplicity  greater  than  unity),  at  which  the  Jacobian  vanishes.  This  is  accomplished  by 
constructing  a  new  system  that  reduces  by  one  the  multiplicity  of  nonsimple  solutions.  Except  in 
degenerate  cases,  this  approach  can  be  applied  sequentially  until  no  roots  have  multiplicity 
exceeding  one. 

The  algorithm  discussed  above  (and  in  more  detail  in  a  preprint  included  in  this  report)  has 
been  used  in  determining  the  coefficients  of  low-dimensional  ordinary  differential  equations 
systems  describing  scalar  transport  in  a  turbulent  plane  mixing  layer  (see  §2.2). 

2.2.  Low-Order  Models  of  Scalar  Transport  in  a  Turbulent  Plane  Mixing  Layer 

In  joint  work  with  Ari  Glezer  and  Martin  Brooke  at  Georgia  Tech,  we  have  demonstrated  the 
utility  of  the  proper  orthogonal  decomposition  (POD)  as  a  data  compression  technique  for  efficient 
representation  of  massive  amounts  of  spatio-temporal  passive  scalar  data  acquired  in  a  turbulent 
plane  mixing  layer.  We  show,  by  comparison  of  computation  to  experiment,  that  such 
compressed-data  representations  are  useful  in  "multi-step-ahead"  prediction  of  scalar  transport  in  a 
turbulent  plane  mixing  layer.  Such  capability  will  allow  for  drastic  reductions  in  the  dimensionality 
of  the  control  scheme. 

In  our  case,  the  application  of  ultimate  interest  is  phase  correction  in  aero-optic  flows  in  shear 
flows,  in  which  large  coherent  structures  significantly  influence  temperature  and  density,  and 
hence  index  of  refraction  distributions.  A  question  arises  as  to  whether  evolution  of  the  coherent 
structures  can  be  manipulated  in  such  a  way  that  the  optical  phase  distortion  can  be  predicted  or 
controlled.  While  entrainment  of  irrotational  fluid  in  turbulent  shear  flows  is  affected  by  large- 
scale  motions,  molecular  mixing  (and  hence  reduction  of  index  of  refraction  gradients)  ultimately 
takes  place  at  the  smallest  scale.  The  traditional  approach  to  control  of  mixing  through 
manipulation  of  global  instability  modes  of  the  base  flow  depends  on  the  classical  cascading 
mechanism  to  transfer  energy  from  the  large  coherent  structures,  whose  evolution  is  being 
manipulated,  to  the  scales  at  which  molecular  mixing  occurs.  Although  mixing  at  the  smallest 
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scales  is  coupled  to  control  of  large  coherent  structures,  more  efficient  control  of  mixing  might  be 
achieved  by  direct  control  of  both  large-  and  small-scale  mixing  processes. 

A  paper,  to  be  submitted  to  the  Journal  of  Fluid  Mechanics,  describes  our  implementation  of 
the  POD  technique  in  representing  massive  amounts  of  spatio-temporal  temperature  data  in  this 
system.  In  that  paper,  we  also  show  how  several  neural  network-based  prediction  methods  are 
able  to  predict  evolution  of  the  temperature  field  based  on  past  data.  We  also  show  how  another 
approach,  based  on  a  low-order  nonlinear  ODE  system  constructed  from  the  data,  attacks  the  same 
problem.  The  implications  of  the  results  for  flow  control  are  discussed  in  a  final  section. 
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Simulation  of  extraction  of  velocity  from  passive  scalar  data 
in  a  two-dimensional  diverging  channel  flow 

Bonnie  N.  Carpenter®*  and  Arne  J.  Pearlstein^’* 

Department  of  Mechanical  and  Industrial  Engineering,  University  of  Illinois  at  Urbana-Champaign, 

1206  Vi/est  Green  Street,  Urbana,  Illinois  6 1 80 1 

(Received  6  June  1995;  accepted  30  May  1996) 

A  technique  for  determining  /? -dimensional  {n  =  2,3)  velocity  fields  from  measurement  of  n  —  1 
passive  or  reactive  scalars  has  recently  been  developed  [Phys,  Fluids  7,  754  (1995)].  The  method 
utilizes  n  linear,  first-order,  uncoupled  hyperbolic  equations  for  n  velocity  components,  derived 
from  transport  equations  for  «-l  scalars  and  conservation  of  mass.  The  velocity  is  determined  by 
integration  along  the  characteristic  curves  defined  by  the  hyperbolic  equations.  Here  we  present  the 
results  of  a  computational  proof-of-concept  study  for  a  steady,  two-dimensional,  diverging  channel 
flow.  We  consider  the  effects  of  the  resolution  of  the  grid  (and  the  shape  of  its  elements)  on  which 
the  scalar  is  known  and  its  derivatives  are  approximated,  the  integration  step  size  used  to  calculate 
the  velocity  components  along  characteristic  curves,  and  the  effects  of  multiplicative  noise 
introduced  into  the  scalar  field.  The  results  show  that  extraction  of  the  velocity  by  integration  of  the 
equations  along  characteristics  is  stable,  and  that  the  techniques  proposed  for  removal  of 
singularities  are  effective.  For  steady  flows,  we  show  that  noise  in  the  scalar  field  measurements  can 
be  dealt  with  by  extracting  the  velocity  from  the  mean  of  several  noisy  scalar  fields.  ©  1996 
American  Institute  of  Physics.  [51070-6631(96)03109-1] 


I.  INTRODUCTION 

Several  techniques  for  determining  velocity  fields  from 
measurements  of  a  single  scalar  have  been  proposed.^"^  In  a 
recent  paper,^  we  have  discussed  the  attributes  of  these  ap¬ 
proaches  and  presented  a  different  method,  requiring  mea¬ 
surement  of  - 1  passive  or  reactive  scalars  to  uniquely  de¬ 
termine  an  n -dimensional  {n-2  or  3)  solenoidal 
(divergence-free)  velocity  field. 

As  discussed  in  Ref.  5,  our  method  differs  from  previous 
approaches,  in  that  it  recovers  the  exact  velocity  field  in  the 
limit  of  complete  data  (i.e.,  perfect  spatial  and  temporal  reso¬ 
lution  of  noise-free  scalar  fields).  Since  experimental  data  is 
always  acquired  with  finite  spatial  and  temporal  resolution, 
and  is  always  corrupted  to  some  degree  by  noise,  it  is  impor¬ 
tant  to  understand  how  our  method  is  affected  by  these  prop¬ 
erties  of  real  data,  and  to  develop  techniques  for  minimizing 
error  in  the  extracted  velocity  field. 

In  this  paper,  we  address  the  issues  of  spatial  resolution 
and  noise  in  the  context  of  a  steady  two-dimensional  diverg¬ 
ing  channel  flow.  Flow  in  a  diverging  channel  is  an  excellent 
testbed  for  developing  an  understanding  of  noise  and  resolu¬ 
tion  effects,  since  the  geometry  admits  an  exact  solution  of 
the  Navier-Stokes  equations.^  Moreover,  the  divergent  na¬ 
ture  of  the  flow  is  reminiscent  of  a  plane  mixing  layer,  a 
generic  flow  of  interest  in  a  number  of  applications.  Finally, 
for  simple  boundary  conditions,  this  flow  gives  rise  to  a  sin¬ 
gularity  in  the  hyperbolic  equations  along  a  curve,  thereby 
allowing  evaluation  of  the  technique  proposed  for  integration 
tlvough  such  a  singularity.^ 

In  Sec.  II,  the  analysis  and  working  equations  derived 


•^Present  address:  Aerospace  Corporation,  El  Scgundo,  California  90245. 
‘’^Corresponding  author:  Telephone:  (217)  333-3658;  fax:  (217)  244-6534; 
Electronic  mail:  ame@ajpiris.mc.uiuc.edu 


earlier^  for  extraction  of  a  two-dimensional,  solenoidal  ve¬ 
locity  field  from  a  single  scalar  field  are  briefly  reviewed.  In 
Sec.  Ill,  we  discuss  the  diverging  channel  flow  and  suitable 
scalar  boundary  conditions.  In  Sec.  IV,  the  procedure  for 
extracting  the  velocity  field  from  the  resulting  computed  sca¬ 
lar  field  is  described.  The  results  of  the  extraction  are  pre¬ 
sented  in  Sec.  V;  discussion  and  conclusions  follow  in  Sec. 
VI. 


II.  DETERMINATION  OF  TWO-DIMENSIONAL 
SOLENOIDAL  VELOCITY  FIELDS  FROM 
MEASUREMENTS  OF  ONE  SCALAR 

In  this  section,  we  present  the  working  equations  derived 
earlier  ^  for  extraction  of  a  two-dimensional,  solenoidal  ve¬ 
locity  field  from  measurements  of  a  single  scalar.  Our  tech¬ 
nique  is  based  on  recognizing  that  the  transport  equation, 

dS 

—  +  u-VS=aV^S  +  R,  (1) 

ot 

for  a  scalar  S  with  diffusivity  a  is  linear  in  the  velocity 
components,  and  that  it  and  the  solenoidal  continuity  equa¬ 
tion, 

V.u=0,  .  (2) 

constitute  two  linear  equations  in  two  velocity  components. 

The  formulas  will  be  displayed  in  a  Cartesian  coordinate 
system.  (Results  for  general  curvilinear  systems  are  shown  in 
Ref.  5.)  The  second  term  on  the  right-hand  side  of  (1)  allows 
lor  the  scalar  to  react,  with  kinetics  that  are  essentially  arbi¬ 
trary  as  long  as  they  are  known  and  the  rate  of  reaction 
depends  solely  on  S,  Xj  y,  and  t.  If  temperature  is  the  scalar 
for  which  full-field  measurements  are  available,  then  (1)  will 
be  replaced  by  an  energy  equation  with  a  replaced  by  the 
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thermal  diffusivity  k,  and  the  reaction  rate  set  to  zero.  In  the 
nonreactive  case,  in  two  dimensions,  (1)  can  be  written  as 


dS  dS  dS  . 

dt  dx  dy 


(3) 


where  is  the  two-dimensional  Cartesian  Laplacian,  and  u 
and  V  are  the  x  and  y  components  of  u.  Equation  (3),  with 
the  two-dimensional  form 


du  dv 
dx  dy 


=  0, 


(4) 


of  (2)  constitutes  a  pair  of  linear  equations  in  u{x,y,t)  and 
v(x,y,l).  Solving  (3)  for  v  and  substituting  into  (4)  yields 


VSxV 


f  G  ] 

dy\dSldy] 


(5a) 


while  solution  of  (3)  for  u  and  substitution  into  (4)  yields 


VSxV 


(— ) 

dS/dx I 


dx  \dS/dx  ' 


(5b) 


where 


,  dS 


(6) 


Equations  (5a)  and  (5b)  constitute  a  pair  of  linear,  first- 
order,  uncoupled  hyperbolic  equations  in  the  velocity  com¬ 
ponents  u  and  V.  Given  the  scalar  field  data,  the  solution  of 
(5a)  and  (5b)  is  made  unique  by  specification  of  appropriate 
boundary  conditions  on  u  and  v.  Equations  (5a)  and  (5b)  can 
then  be  integrated  along  characteristic  curves,  which  are  the 
same  for  both  equations.  No  time  derivatives  of  the  velocity 
components  appear  in  (5a)  and  (5b),  so  that  no  initial  data  on 
the  velocity  field  are  required.  As  this  is  not  an  initial  value 
problem,  errors  in  the  determination  of  the  velocity  field  will 
not  grow  temporally. 

Equations  (5a)  and  (5b)  have  an  apparent  singularity 
when  either  component  of  the  gradient  vanishes.  As  shown 
in  Ref.  5,  these  apparent  singularities  are  “removable’*  in  all 
cases  except  when  V5  vanishes  identically  over  an  area.  The 
singularities  are  removable  in  the  sense  that  the  velocity  field 
can  be  determined  at  points  or  along  curves  where  one  or 
both  components  of  V5  vanish,  and  that  U'VS  =  0  over  an 
area  poses  no  difficulty  unless  V5=0  identically.  If  only  one 
component  of  VS  vanishes,  then  the  integration  can  be  per¬ 
formed  for  one  velocity  component  using  the  nonsingular 
equation  of  (5a)  and  (5b),  and  the  other  velocity  component 
can  be  determined  from  the  scalar  transport  equation  (3).  If 
both  components  of  VS  vanish  at  a  point  or  on  a  curve,  then 
u  can  be  determined  as  follows.  Using  index  notation,  we 
rewrite  (3)  as 


(9S 

Uj 

^  dxj 

Taking  the  gradient  of  both  sides  of  (7),  we  obtain 
d^S  dG 


^  dXi  dXj  dXi' 


1  =  1,2, 


(7) 


(8) 


FIG.  1.  Flow  in  a  diverging  channel.  Velocity  field  shown  schematically. 


when  VS  vanishes.  At  each  point,  (8)  is  a  system  of  two 
linearly  independent,  nonhomogeneous  algebraic  equations 
for  the  velocity  components. 


III.  THE  SCALAR  FIELD 

To  better  understand  the  properties  of  the  proposed  pro¬ 
cess  for  determining  the  velocity  field  from  passive  scalar 
measurements,  a  computational  proof-of-concept  study  has 
been  performed  using  advection  of  a  passive  scalar  in  steady 
two-dimensional  diverging  channel  flow.  We  can  compare 
the  velocity  extracted  from  the  computed  scalar  field  to  the 
Jeffery-Hamel  solution,^  an  exact  solution  of  the  Navier- 
Stokes  equations.  This  involves  three  steps:  calculation  of 
the  Jeffery-Hamel  flow,  calculation  of  the  scalar  field  using 
the  Jeffery-Hamel  solution,  a  scalar  transport  equation,  and 
appropriate  boundary  conditions,  and  extraction  of  the  veloc¬ 
ity  field  from  the  computed  scalar  field  by  integration  along 
characteristics. 

One  goal  of  this  work  is  to  investigate  the  effects  of 
spatial  resolution  and  noise  on  the  accuracy  of  the  extracted 
velocity  components.  In  this  steady  two-dimensional  case, 
the  sources  of  error  in  the  extracted  velocity  field  include  the 
spatial  resolution  of  the  scalar  field  used  to  approximate  de¬ 
rivatives  of  the  scalar  in  (5a)  and  (5b),  and  the  integration 
step  size  along  the  characteristics  used  in  extracting  the  ve¬ 
locity  components.  We  also  consider  the  effect  of  multipli¬ 
cative  noise  (e.g.,  shot  noise  in  full-field  optical  measure¬ 
ment  of  temperature  by  laser-induced  fluorescence,  etc.) 
corrupting  measurements  of  the  scalar  field. 

A.  Exact  velocity  field 

We  consider  steady  two-dimensional  flow  in  a  diverging 
channel,  as  shown  in  Fig.  1.  In  this  case,  there  is  an  exact 
solution  oY  the  Navier-Stokes  equations,  referred  to  as 
Jeffery-Hamel  flow,^  with  the  velocity  field  depending  on 
the  half  angle  <t>^  and  the  Reynolds  number,  Rt=UQr/v, 
where  UQ(r)  is  the  velocity  along  the  central  streamline  {<f> 
=0).  In  cylindrical  coordinates,  this  purely  radial  flow  can  be 
written  as 
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yF(<}>) 

U=  — - —  e,, 

where  F  is  an  even  function  of  (f>  given  by 


(9) 


F(<^)  =  Re-6 

with  k  found  from 

i+k^  2 

3it^(l+2/Re)“*” 

for  Re<Re*((^o),  or  by 


l  +  Re/2\ 

”TTP 


\k^  sn^ 


l+Re/2 

"T+F 


1/2 


4>,k 


l+Re/2\ 

T+F 


<^o.* 


F{<t>)  = 


Re-  q  +  (Re+  q)cr\[{2q/3y'^<f>,k] 


1  +cn[(2(y/3)'^<^,/:] 
with  k  and  q  defined  by 


6  +  Re(5-4Jt^) 


6  +  Rt{\+4k^) 


2  +  Re\''^ 


2k^-\ 


and 


3(1+ Re/2) 
2k^-  \  ~  ' 


(10a) 


(10b) 


(11a) 


(11b) 


(11c) 


characteristics  originating  on  the  wall  emanate  into  the  flow 
or  in  the  opposite  direction  out  of  the  channel.  If  a  charac¬ 
teristic  originating  on  the  wall  does  not  emanate  into  the 
flow,  then  the  direction  of  integration  along  the  characteristic 
must  be  reversed.  This  can  be  accomplished  easily  by  a 
change  of  variables.  Letting  u '  =  —  «  and  i;  ^  —  u  and  sub¬ 

stituting  into  (12a)  and  (12b)  yields 


0  i^^S)Syy 


^y^xx\  , 


(^^5)5, 


(13a) 


(13b) 


The  net  result  is  to  change  the  signs  of  the  right-hand  sides, 
so  that  characteristics  that  did  not  previously  emanate  into 
the  flow  will  now  do  so. 

A  simple  choice  for  the  boundary  conditions  on  the  sca¬ 
lar  is 


for  Re>Re*(<;^^o),  where  sn  and  cn  are  Jacobian  elliptic  func¬ 
tions,  and  Re*(5®)  =  684.^  Bisection  is  used  to  determine  k  in 
(10b)  or  (11b).  Once  k  is  known,  the  velocity  field  can  be 
found  using  Eq.  (10a)  for  Re>Re*(</>o)  or  Eqs.  (lla)-(llc) 
for  Re<Re*(0j. 


B.  Computed  scalar  field 


We  use  the  exact  velocity  field  to  calculate  the  scalar 
field.  Before  choosing  boundary  conditions  for  the  scalar,  it 
is  useful  to  consider  the  nature  of  the  characteristic  equations 
for  steady  two-dimensional  solenoidal  flow.  In  recognition  of 
the  fact  that  experimental  data  will  frequently  be  available  on 
domains  whose  boundaries  do  not  coincide  with  constant 
coordinate  curves  or  surfaces  of  simple  coordinate  systems, 
we  will  extract  the  velocity  components  using  a  Cartesian 
system,  rather  than  the  radial  coordinates  that  are  natural  for 
the  geometry  at  hand.  Thus,  (5a)  and  (5b)  reduce  to 


=  a 


(V^5)5. 


=  -a  (V25),-- 


5. 


(12a) 


(12b) 


where  subscripts  denote  partial  derivatives.  From  (12a)  and 
(12b)  one  can  see  that  the  characteristic  curves  depend  solely 
on  the  .X  and  y  derivatives  of  the  scalar.  Therefore,  in  order 
for  the  characteristics  to  propagate  away  from  a  wall  (at 
<f>-  ±  (/>^ ,  on  which  «  =  u  =  0),  the  scalar  on  at  least  one  wall 
must  be  nonuniform.  Note  that  if  and  Sy  are  nonzero, 


( 


5(r,-(^,)=|  5,  +  (^^|(S2-Si), 


52,  r2=^r, 


5(r,</>,)  =  5i. 


ri^r^r2. 


(14a) 

(14b) 


However,  due  to  the  discontinuity  in  the  radial  derivative  of 
5  on  the  (/>=  -  <t>o  wall,  the  corresponding  scalar  field  is  not 
sufficiently  differentiable  to  allow  computation  of  the  third 
derivatives  needed  in  (12a)  and  (12b).  To  obtain  a  suffi¬ 
ciently  smooth  scalar  field,  we  will  require  that  Dirichlet 
conditions  on  S  be  continuously  differentiable.  Such  a  set  of 
boundary  conditions  is 


S(r,-(^J  =  5i  +  A5 


tanh  i7[(r-ri)/rj]  +  tanh  b 
1  +  tanh  b 


=  5^/(0, 


5(r,<^J  =  5i, 


(15a) 

(15b) 


for  rj>0,  fc>0. 

At  this  point,  we  nondimensionalize  by  defining 
^=(5-5i)/AS,  and  4>^)I2<I)^.  The  dimensionless 

scalar  transport  equation  and  boundary  conditions  are 


F{(f>)de_  \  d  /  1 


e(e,0)= 


tanh  b{^—  1  )  +  tanh  b- 
1  +  tanh  b 


(16) 

(17a) 


0(M)  =  O,  (17b) 

where  F{<{>)-F{2(}>q<I)-  and  a=vla  is  the  ratio  of  mo¬ 
mentum  and  scalar  diffusivities.  We  see  that  vanishes 
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at  ^=0,  and  approaches  unity  as  The  parameter  b 

characterizes  the  sharpness  of  the  transition  of  the  scalar  dis¬ 
tribution  on  the  wall  at  We  represent  the  solution  of 
d{l^)  as 

oc 

n  -  1 

where  g  will  be  chosen  to  satisfy  the  boundary  conditions 
(17a)  and  (17b),  and  each  term  in  the  summation  will  satisfy 
homogeneoi^  boundary  conditions  on  the  walls. 

If  g(^y<f))  and  the  Fourier  expansion  each  satisfy  the 
scalar  transport  equation  at  the  walls,  the  scalar  distribution  6 
will  be  sufficiently  differentiable  there.  Since  F=0  at  the 
walls,  we  require 


Setting 

one  can  show  that 

x[2fef  tanhfc(^-l)-l], 


g3(^)=k2(^).  (21c) 

g,(e)=-M^)-fg2(^)-  (21d) 

In  light  of  (18),  (16)  becomes 

“  [  Ii-F(^)<7\  /  «Tr\M  1  ,  _ 

F{^)dg  1  .?  Idg]  1  Si 
^  di  44>ie  dcP^' 

where  primes  deno^  differentiation  withjrespect  to  Multi¬ 
plying  by  sinm7r<^,  integrating  from  4>-0  to  0=1,  rear¬ 
ranging  terms,  and  approximating  the  radial  derivatives  of 
each  function  by  second-order  central  differences,  we  ob¬ 
tain 


(19) 

(20) 

m=l 

/  a  \ 

,  (23) 

(21a) 

=  2(A^)2^-F„(f,.)+-C„(e,)J 

where  +  «n(^y)  is  approximated  by  and 

Amn-l  F(0)sin /i7r0  sin  m7r0  f/0,  (24a) 

Jo 


fifl  ^  /  o»g\  1  d^g]  ^  -  (-1)^ 

Jo  1?  ^  I  ^  ^  ^ 


_  _  (-l)"-^'  +  l  fl  d  j  ^dgo 


(-1)"^'  [1  d  I  ^dg^\]  /2[(-l)"-l3  (-l)"^'\  1  d  /  ^2] 

077  nV^  ’’’  riTT 

1  /(-1)"^'  +  1\  /6(-l)'’  (-1)”^‘\[1  d  I  dg:,\]  3  (-I)"-"' 

^ 2^ (  '.V- « 31  (««))  / 2^ 


C„(0=  f  17  F(<p)sin  mr<p  d^ 
Jo  <?€  , 


+  — ^  f  ^F(<p)sin  mT(f>  d<p 

d^\^h 


0^F{0)sin  /2  7r0  d<f>. 


F(0)sin  niT<f>  d<(> 


^■“rr  0F(0)sin  n7r0 
L  Jo 


with  a^n  =  ai  =  0.  The  infinite  series  in  (23)  is  then  truncated 
to  N  terms.  We  choose  to  be  large  enough  so  that  the 
scalar  field  at  ^  is  approximately  linear  in  <f>. 

This  system  of  linear  algebraic  equations  for  the  coeffi¬ 
cients  a{  was  solved  using  standard  numerical  linear  algebra 
software.  The  scalar  distribution  ^  at  a  given  point  is  deter¬ 
mined  from  (18)  using  a  finite  number  of  terms  and  interpo¬ 
lation  between  the  discrete  values  of  ^  at  which  the  coeffi¬ 
cient  functions  have  been  approximated. 
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IV.  EXTRACTION  OF  VELOCITY  FIELD  FROM  THE 
COMPUTED  SCALAR  FIELD 


The  final  step  is  to  calculate  the  velocity  by  the  method 
of  characteristics.  The  nondimensional  forms  of  (12a)  and 
(12b)  are 

0re^.\ 


1 

Pe 


-  ^  1 

u 

{^^e)ey 

z] 

By 

)’ 

Bye„^ 

ji; 

'■  B,  t 

) 

(25a) 


(25b) 


where  Pe  =  c7-Re  is  the  Peclet  number.  To  integrate  these 
along  a  characteristic  we  rewrite  (25a)  and  (25b)  as 


— 

ds 


(26a) 


(26b) 


du 


e,e„\ ,  1 


(9v  Pe 


+  s: 


iv^e)eyy 

e. 


dv 


6,6, 


e,  Pe 


{V^e)e, 

0> 


(26c) 


(26d) 


Equations  {26a)-(26d)  are  integrated  along  characteristics  by 
a  Runge-Kutia  method  with  constant  step  size.  The  x  and  y 
derivatives  can  be  obtained  from  ^  and  <f>  derivatives  by  the 
chain  rule. 

Recalling  that 
v 

S  a„(^)s\n  nir$+g{^,^),  (27) 


we  locally  approximate  the  coefficient  functions 
their  first  three  radial  derivatives  by  fourth-degree  Lagrange 
polynomials  fitting  the  values  a{.  To  illustrate  this,  consider 
integration  along  a  characteristic  through  a  point  {x^yy^),  as 
shown  in  Fig.  2.  The  radial  grid  line  closest  to  is  7*3, 

and  a  Lagrange  polynomial  is  fitted  through  the  five  radial 
values  U  =  7i  J2  J3  J4  J5)  approximate  a„(0  at 
Radial  derivatives  of  a„(^  are  approximated  by  de¬ 
rivatives  of  the  polynomial.  Computation  of  derivatives  by 
this  method  provides  a  good  test  of  the  efficacy  of  “singu¬ 
larity  removal”  by  the  procedure  discussed  at  the  end  of  Sec, 
II.  However,  this  method  for  obtaining  the  necessary  deriva¬ 
tives  of  the  scalar  provides  more  accurate  information  than 
one  would  expect  in  experiment,  since  the  azimuthal  deriva¬ 
tives  are  computed  from  an  analytical  expression  rather  than 
by  numerical  differentiation  of  spatially  discrete  data.  For 
this  reason,  the  scalar  was  also  evaluated  at  locations  on  a 
rectangular  grid  using  (27),  as  discussed  in  Sec.  V  A.  The 


FIG.  2.  Relationship  of  radial  grid  lines  to  a  characteristic  curve. 


derivatives  of  6  found  from  the  chain  rule  were  then  com¬ 
puted  by  second-order  accurate  central  difference  approxi¬ 
mations,  or  by  second-order  accurate  lopsided  or  one-sided 
difference  approximations  where  necessary. 

The  integration  is  initiated  at  a  number  of  points  along 
the  wall  and  continued  until  the  characteristic 

curves  leave  the  computational  domain  at  ^ .  This  provides 
a  velocity  field  at  points  on  the  characteristics  (cf.  Fig.  3). 
Velocity  profiles  at  a  given  jc  position  can  be  obtained  from 
the  value  on  each  characteristic  crossing  a  constant-jr  line. 

As  discussed  in  Sec.  II  and  in  more  detail  in  Ref.  5,  the 
vanishing  of  either  component  of  leads  to  an  apparent 
singularity.  In  the  present  case,  this  corresponds  to  the  van¬ 
ishing  of  dO/dx  {dd/dy  is  nonzero  throughout  the  domain). 
We  have  examined  four  ways  to  extract  u  and  v . 

(1)  Extract  u  and  i;  from  (25a)  and  (25b),  except  when 
use  of  a  dimensionless  scalar  transport  equation  analogous  to 
(3)  is  necessary  to  pass  over  an  apparent  singularity  in  (25b), 
as  discussed  in  Sec.  II. 

(2)  Extract  u  from  (25a)  and  compute  v  from  the  analog 
of  (3). 

(3)  Extract  v  from  (25b)  and  compute  u  from  the  analog 
of  (3),  except  at  points  where  d6/dx  vanishes,  where  the  ex¬ 
traction  proceeds  according  to  method  (2), 

(4)  Extract  one  velocity  component  using  the  hyperbolic 
equation  (25a)  or  (25b)  in  which  the  denominator  of  the 
second  term  in  the  second  factor  on  the  right-hand  side  has 
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the  largest  magnitude.  The  other  velocity  component  is  ob¬ 
tained  from  the  analog  of  (3).  As  discussed  in  Sec.  V  A,  we 
expect  the  fourth  approach  to  extract  the  most  accurate  ve¬ 
locity  field. 

Several  noise  processes  can  corrupt  optical  measure¬ 
ments  of  the  scalar,  and  can  be  expected  to  lead  to  errors  in 
the  extracted  velocity  field.  Here,  we  consider  multiplicative 
noise  corrupting  M  measurements  of  a  steady  scalar  field  on 
a  rectangular  grid,  so  that  each  measured  field  is  related  to 
the  actual  field  by 

^deiLied^^i  7)  “  ^actua](‘^j  ^yj)[  ^ 

(l^m^M),  (28) 

where  77  is  the  maximum  amplitude  of  the  noise  and  Z>  is  a 
random  number  with  \D\^\  at  each  point  on  the  grid.  If  the 
temperature  field  were  computed  M  times  and  then  averaged, 
we  would  obtain 

{  ^deiccied(^J  ^actual! -^i 

+  ,yj)  e^^(Xi  ,yj)),  (29) 

where  {  )  indicates  an  ensemble  average.  For  large  Af ,  (29) 
reduces  to 

(  ^dciecied(*^i  ^actual(-^/  yyj)y  (^^) 

since  (D)— ^0  for  a  sufficiently  large  number  of  realizations. 
An  alternative  approach  would  be  to  extract  the  velocity 
from  each  of  the  M  computed  scalar  fields  and  then  average 
the  extracted  velocities.  Besides  being  more  cumbersome, 
this  method  typically  yields  the  wrong  result  (even  in  the 
limit  for  nonzero  77,  since  the  noise  appears  nonlin- 

early  in  the  hyperbolic  equations  for  u  and  u. 


V.  RESULTS 

Before  considering  the  extracted  velocity  fields,  we  first 
show  that  the  scalar  field  computation  converges  as  J  and  N 
increase.  For  all  results  shown  here,  Re=2,  cr=7  (corre¬ 
sponding  to  the  Prandtl  number  for  water),  and  b-2.  Figure 
4  shows  the  absolute  difference  on  the  line  x~2  between  6 
calculated  for  160  and  N=48  and  6  calculated  at  lower 
resolutions.  The  calculation  of  the  scalar  (hereinafter  referred 
to  as  temp>erature)  distribution  converges  rapidly,  and  the 
absolute  difference  between  ^7=i60.n=48  ^7*=40jv=i2 

less  than  10”^.  Therefore,  we  will  consider  the  “exact”  tem¬ 
perature  field  to  be  converged  when  computed  with  40  radial 
grid  points  and  12  azimuthal  modes. 

Most  present  optical  detectors  utilize  identical  rectangu¬ 
lar  (usually  square)  elements  in  a  regular  array,  with  implicit 
signal  averaging  over  the  area  of  each  element.  Thus,  the 
temperature  data  will  normally  present  themselves  on  a  rect¬ 
angular  grid,  uniformly  spaced  in  both  Cartesian  directions. 
Since  in  our  simulation  computation  of  6  is  not  performed  on 
a  rectangular  grid,  interpolation  is  required  to  effect  a  rect¬ 
angular  presentation  of  the  data.  This  constitutes  an  addi¬ 
tional  source  of  error  in  our  simulated  extraction  of  u. 

Extraction  of  u  from  data  on  a  rectangular  grid  requires 
approximation  (e.g.,  by  finite  differences)  of  seven  x  and  y 
derivatives  of  d.  The  magnitude  of  the  error  introduced  into 
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FIG.  4.  Absolute  differences  at  jr  =  2.0  between  the  scalar  field  cornpuied 
with  J=160,  N-48,  and  those  computed  with  less  accuracy,  (a)  {J,N) 
=  (10.3),  (20,6),  (30,9),  (40,12).  (b)  (J.N)  =  (20.6).  (30,9),  (40,12). 
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u  will  depend  on  the  noise  level,  grid  size,  the  difference 
schemes  used  to  approximate  the  derivatives,  and  the  algo¬ 
rithm  and  step  size  used  to  integrate  the  hyperbolic  equations 
along  the  characteristics. 

A.  Noise-free  results 

We  first  illustrate  the  extraction  of  u  and  the  nature  of 
the  apparent  singularity  and  its  removal,  in  the  absence  of 
noise  in  0  (77=0). 

We  begin  by  extracting  u  using  azimuthal  and  radial 
derivatives  of  6  and  the  chain  rule  to  approximate  x  jmd  y 
derivatives.  We  use  (27)  to  evaluate  the  azimuthal  {<f>)  de¬ 
rivatives  of  6,  and  fourth-degree  Lagrange  polynomials  to 
evaluate  radial  derivatives  at  points  on  the  characteristics. 
This  allows  us  to  extract  u  with  less  error  in  approximating 
the  derivatives  of  6  than  would  be  incurred  if  interpolated 
values  of  0  were  numerically  differentiated  on  a  rectangular 
grid. 

Integrating  the  hyperbolic  equations  for  u  and  v  without 
removing  the  singularity  by  the  methods  discussed  in  Sec.  II 
and  Sec.  IV,  we  obtain  the  results  shown  in  Fig.  5.  The 
integration  step  size  is  A5  =  10“^.  Points  in  Figs.  5(a)  and 
5(b)  represent  values  of  the  x-  and  y-velocity  components, 
respectively,  on  characteristics  crossing  the  line  x=1.8.  Fig- 
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FIG.  5.  Here  ;r=  1.8,  Re=2,  cr^l,  i>  =  2,  7=40,  N=12,  and  A5=10”^  Here  u  and  v  exiracied  directly  from  hyperbolic  equations.  Filled  symbols  denote  the 
boundary  data,  (a)  Extracted  (O)  and  exact  (— )  jt  component;  (b)  y  component  of  velocity;  (c)  relative  error  in  the  jc  component;  (d)  y  component  of  velocity. 


ures  5(c)  and  5(d)  the  relative  errors  in  u  and  v,  respectively. 
At  x=1.8,  the  relative  error  in  v  is  much  larger  near  the 
nonisothermal  wall  than  near  the  isothermal  one.  The  char¬ 
acteristic  curves  on  which  the  error  in  v  is  large  are  those 
that  have  crossed  the  curve  SO/dx—0  shown  in  Fig.  6.  (Note 
that  d6ldy<0  everywhere  in  the  domain.)  This  is  what  is 
expected  for  extraction  without  removal  of  the  singularity  in 
(25b).  By  integrating  the  hyperbolic  equation  for  u  and  com¬ 
puting  V  from  the  energy  equation,  one  obtains  the  results 
shown  in  Fig.  7,  in  which  the  relative  errors  in  v  are  much 
lower  than  those  shown  in  Fig.  5. 

The  relative  error  for  u  in  Fig.  7(c)  increases  to  approxi¬ 
mately  5%  on  the  characteristics  closest  to  the  isothermal 
wall.  The  values  of  i;  on  these  characteristics,  computed 
from  u  and  the  energy  equation,  are  in  error  by  as  much  as 
10%.  To  understand  the  source  of  this  error,  consider  again 
the  characteristic  curves  shown  in  Fig.  3.  The  characteristics 
lying  close  to  the  isothermal  wall  must  traverse  a  greater 
distance  to  reach  x  =  1.8,  but  the  integration  is  stable  and 
errors  do  not  grow  rapidly  along  the  integration  path.  We 
believe  that  the  larger  errors  are  due  to  errors  in  approximat¬ 
ing  radial  derivatives  of  temperature  in  the  upstream  region. 
Finally,  we  note  that  the  relative  error  in  i;  is  approximately 
twice  that  in  u.  This  is  not  surprising  since  v  depends  on  a 
previous  computation  of  u. 

To  more  faithfully  mimic  extraction  of  u  from  experi¬ 
mental  data,  the  lemj)erature  was  computed  at  locations  on  a 
rectangular  grid  using  (27),  and  its  derivatives  were  approxi¬ 
mated  by  finite  differences.  The  resolution  of  the  grid  is  de¬ 
noted  by  (Ig  Jg),  where  Ig  and  Jg  are  the  number  of  x  and  y 
grid  points,  respectively.  With  /g=215  and  7^=210,  we  first 
investigate  the  effect  of  integration  step  size.  Figures  8  and  9 
show  results  for  A5  =  10"*  and  10"^,  respectively.  These  re¬ 


sults  were  obtained  by  integrating  the  hyperbolic  equation 
for  u  and  computing  v  from  the  energy  equation.  The  rela¬ 
tive  errors  in  both  extracted  velocity  components  are  much 
smaller  for  A5  =  10~^  than  for  A5  =  10  ^  For  A5  =  10  ^  (not 
shown),  the  relative  error  decreases  slightly  at  a  few  points 
near  the  isothermal  wall  but  is  otherwise  unchanged.  De¬ 
creasing  £^s  io  lO”"*  gives  no  further  reduction  in  error. 
Therefore  the  remaining  error  is  due  to  approximating  the 
derivatives  of  0  by  finite  differences  using  6  values  interpo¬ 
lated  onto  the  points  of  a  rectangular  grid.  Figures  7  and  9 
show  that  the  maximum  error  in  u  doubles  from  5%  to  10% 
and  that  for  v  doubles  from  10%  to  20%  as  a  result  of  finite- 
difference  approximation  of  the  derivatives.  Thus,  the  errors 
incurred  by  the  integration  (for  A5=10“^)  are  much  smaller 
than  those  due  to  approximation  of  temperature  derivatives 
on  a  grid.  Since  the  relative  errors  in  the  velocity  compo- 


FIG.  6.  Sign  of  6^  for  Jeffery-Hamel  flow  (Re=2,  <r=l,  b-2). 
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FiG.  7.  Here  1.8,  Re  =  2,  tT=7,  =  J=40,  N=12,  and  ^s=\0~^.  Here,  u  is  extracted  directly  from  the  hyperbolic  equation;  v  is  computed  from  u  and 

the  scalar  transport  equation.  Filled  symbols  denote  the  boundary  data,  (a)  Extracted  (O)  and  exact  ( — )  x  component;  (b)  y  component  of  velocity;  (c)  relative 
error  in  the  x  component;  (d)  y  component  of  velocity. 


nenls  are  virtually  the  same  for  A5  — 10  10  and  10 

the  remainder  of  the  calculations  will  be  performed  with 
A5=10"'1 

We  next  investigate  the  effects  of  spatial  resolution  and 
element  shape  on  extraction  of  the  velocity  components.  Fix- 
ing  A5-10'^  7  =  40,  and  A^=12  as  discussed  above,  and 
considering  grids  with  (/^  ,7^)  =  (200,34),  (300,52),  (400,68), 


and  (5(X),86),  we  obtain  the  results  shown  in  Fig.  10.  Again, 
these  results  are  obtained  by  integrating  the  hyperbolic  equa¬ 
tion  for  u  and  computing  i;  from  the  energy  equation.  Fig¬ 
ures  10(a)(i,ii)  show  that  the  errors,  especially  near  the  iso¬ 
thermal  wall,  are  very  large  for  this  coarse  grid.  When  the 
resolution  is  increased  to  (300,52),  the  error  decreases  sig¬ 
nificantly  [Figs.  10(b)(i,ii)]  but  is  still  quite  large  (>20%)  at 


FIG.  8.  Here  x=1.8,  Re=2.  <r=7,  J=40.  A^=12,  /,=215,  and  7^=210.  Here,  u  is  extracted  direcUy  from  the  hyperbolic  equation;  i;  is 

computed  from  u  and  scalar  transport  equation.  Filled  symbols  denote  the  boundary  data,  (a)  Extracted  (O)  and  exact  {  )  x  component;  (b)  y  component  of 
velocity;  (c)  relative  error  in  the  x  component;  (d)  y  component  of  velocity. 
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nC.  9.  Here  jr=1.8,  Re  =  2,  a-=7»  b  =  2,  7=40.  N=12,  AJ  =  10■^  /,  =  215,  and  7^  =  210.  Here,  u  is  extracted  directly  from  the  hyperbolic  equation;  v  is 
compuicd  from  u  and  the  scalar  transport  equation.  Filled  symbols  denote  the  boundary  data,  (a)  Extracted  (O)  and  exact  ( — )  x  component;  (b)  y  component 
of  velocity;  (c)  relative  error  in  the  x  component;  (d)  y  component  of  velocity. 


some  points.  When  the  resolution  is  increased  further  to 
(400,68)  [Figs.  10(c)(i,ii)],  the  error  is  further  reduced  and 
exceeds  20%  only  near  the  isothermal  wall.  For  (5(X),86) 
[Figs.  10{d)(i,ii)],  the  relative  error  in  u  is  always  less  than 
12%,  while  the  error  in  v  exceeds  10%  only  on  the  two 
characteristics  closest  to  the  isothermal  wall.  This  is  due  to 
errors  incurred  in  approximating  temperature  derivatives  at 
and  near  the  upstream  origin  of  these  characteristics.  A  lo¬ 
cally  refined  grid  could  be  used  to  reduce  this  source  of  error. 

Comparison  of  Fig.  9  to  Figs.  10(c)  and  10(d)  shows  that 
the  grid  of  nonsquare  elements  (Fig.  9)  gives  results  with 
accuracy  comparable  to  that  of  a  grid  of  square  elements 
[Fig.  10(c)]  having  about  40%  fewer  points  (45  150  versus 
27  200),  and  that  the  accuracy  is  inferior  to  that  of  a  grid 
with  a  comparable  number  (43  000)  of  square  elements  [Fig. 
10(d)]. 

Four  approaches  for  extracting  the  velocity  components 
from  the  hyperbolic  equations,  and  possibly  the  energy  equa¬ 
tion,  were  discussed  in  Sec.  IV.  The  numerical  results  de¬ 
scribed  above  were  found  with  Approach  2.  Setting  J=40, 
N=12,  and  A5=10‘^  as  before  and  fixing 
(7^  ,7^)  =  (400,68),  we  compare  the  accuracy  of  the  four  ap¬ 
proaches.  Approach  1  extracts  u  and  v  from  the  hyperbolic 
equations  except  when  the  energy  equation  is  necessary  to 
bypass  a  singularity.  The  x  component  is  identical  to  that 
shown  in  Fig.  10(c)(i),  and  the  y  component  is  shown  in  Fig. 
11(b).  As  discussed  above,  the  results  from  Approach  2  are 
'  shown  in  Figs.  10(c)(i,ii).  Those  from  Approach  3  (u  ex--“ 
tracted  from  the  hyperbolic  equation  and  u  computed  from 
the  energy  equation)  are  shown  in  Figs.  11(a)  and  11(b). 
Results  from  Approach  4  [one  velocity  component  extracted 
from  the  hyperbolic  equation  (12a)  or  (12b)  having  the  larger 
denominator  in  the  second  term  of  the  second  factor  on  its 


right-hand  side,  with  the  other  component  computed  from 
the  energy  equation]  are  identical  to  those  shown  for  Ap¬ 
proach  2  in  Figs.  10(c)(i,ii).  This  is  because  \d6/dy\ 
>\d6/6x\  throughout  the  domain,  so  that  Approaches  2  and 
4  are  identical.  Figure  11(b)  shows  that  extraction  of  i;  from 
its  hyperbolic  equation  results  in  large  errors,  since  \d9ldx\  is 
small  [0(10“^)]  throughout  the  domain,  especially  when 
close  to  (but  not  on)  the  curve  \dBldx\=^  shown  in  Fig.  6. 
Figure  11(b)  shows  that  extraction  of  u  from  the  energy 
equation  using  values  of  v  computed  from  its  hyperbolic 
equation  results  in  large  errors  in  u  also.  Therefore  the  most 
accurate  method  for  extracting  the  velocity  is  Approach  4 
(identical  to  Approach  2  for  this  particular  flow).  All  subse¬ 
quent  results  are  obtained  using  Approach  2. 

B.  Results  with  noise 

To  investigate  the  effects  of  noise  on  the  error  in  the 
velocity  extraction,  the  temperature  field  was  corrupted  by 
multiplicative  noise  at  the  grid  points.  The  most  accurate 
temperature  field  available,  7  =  160,  N  =  12,  with  an  integra¬ 
tion  step  size  A5  =  10”^,  was  used  for  all  the  following  re¬ 
sults,  unless  otherwise  specified.  For  noise  amplitudes 
77=10”'^,  10“^  and  10“®,  Figs.  12(a)(i,ii),  12(b)(i,ii),  and 
12(c)(i,ii),  respectively,  show  the  extracted  velocity  for  a 
single  realization  (M  =  l),  with  a  grid  composed  of  non¬ 
square  elements  (7^=215,  7^  =  210).  For  77=10  ^  the  errors 
in  both  velocity  components  are  terge  (>100%).  When  77  is 
reduced  to  10"^,  the  errors  in  both  velocity  components  are 
still  large,  but  significantly  smaller  than  for  77=  10”^.  When 
77  is  reduced  to  10“®,  the  errors  incurred  due  to  noise  are 
insignificant,  and  the  results  are  similar  to  the  noise-ffee 
case. 
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Fir  in  Here  x- 1  8  Re  =  2  <t=7  i>  =  2  J=40.  W=12,  and  Ai=10-2;  extracted  (O)  and  exact  (— )  components,  Here,  u  is  extracted  directly  A* 
TytrS^lic  eqoat;  0  ts  computed  from  «  and  Uie  scalar  transport  equation.  Filled  symbols  denote  the 

cogent,  for  /.  =  200,  J,  =  34;  (b)(i)  x  component;  (b)(ii)  y  component,  for  y,  =  300. 2, =52;  (c)(i)  x  component;  (c)(ii)  y  component,  for  /,-400. 7.  68. 
(d)(i)  X  component;  (d)(ii)  y  component,  for  /,=500,  7,-86. 


If  a  grid  with  square  elements  (7^-400,  68)  is  used 

instead,  the  results  shown  in  Figs.  13(a)  and  13(b)  are  ob¬ 
tained  for  ?;=10'''  and  10■^  respectively.  (In  what  follows, 
both  u  and  i>  were  computed  by  Approach  2,  as  discussed 
above,  and  reference  is  made  to  results  for  both  components, 
although  results  are  shown  only  for  u.)  For  77=10  the 
errors  in  both  velocity  components  are  large,  but  smaller 
than  those  shown  in  Figs.  12(a)(i,ii)  for  nonsquare  elements. 
When  77  is  reduced  to  10"*,  the  errors  in  both  u  and  v  are 
similar  to  those  with  no  noise.  For  77=10”®,  the  results  are 
indistinguishable  from  those  with  no  noise,  and  nearly  indis¬ 


tinguishable  from  those  with  77=10  *.  Comparison  of  the 
noise-free  results  for  nonsquare  and  square  elements  [Figs. 
9(a)  and  10(c)(i),  respectively]  to  the  corresponding  results 
for  77=10"*  [Figs.  12(b)(i)  and  13(b)]  indicates  that  extrac¬ 
tion  of  the  velocity  components  is  less  susceptible  to  noise- 
---related  errors  when  the  elements  of  the  grid  are  square.  It 
was  noted  earlier  that  \deidy\>\d0ldx\  throughout  the  do¬ 
main,  so  that  we  extract  u  from  its  hyperbolic  equation  and 
compute  V  from  the  energy  equation.  Since  00/dy  appears  in 
the  denominators  of  two  terms  in  the  hyperbolic  equation  for 
u,  errors  in  approximating  dOldy  will  be  exacerbated  relative 
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HG.  11.  Here  ;t=1.8,  Re=2,  cr=7.  fc=2.  2=40.  N=12,  Ar  =  10"^  /,=400,  and  2, =68.  Here,  v  is  extracted  directly  from  the  hyperbolic  equation;  u  is 
computed  from  v  and  the  scalar  transport  equation.  Filled  symbols  denote  the  boundary  data,  (a)  Extracted  (O)  and  exact  ( — )  x  component;  (b)  y  component 
of  velocity. 


to  those  for  dO/dx  as  Ay/ Ax  decreases,  since  6  is  corrupted 
by  noise  at  each  grid  point. 

We  can  reduce  the  effect  of  noise  on  the  extracted  ve¬ 
locity  by  averaging,  as  discussed  in  Sec.  IV.  First,  we  con¬ 
sider  extraction  of  u  and  v  from  each  of  M  noisy  tempera¬ 
ture  fields,  followed  by  averaging  the  velocity  fields.  By 
fixing  lg  =  ^00,  ^^=68,  and  the  noise  amplitude 


0 


l7=10“^  and  taking  M  =  l,  100,  we  obtain  the  results  shown 
in  Figs.  14(a)  and  14(b),  respectively.  (For  these  results, 
J=40  and  N=  12.)  The  noise  has  little  effect  on  the  error  in 
u  and  u  for  M  =  l,  or  for  M  =  10  or  30  (not  shown,  but 
indistinguishable  from  the  M  =  l  case).  However,  for 
M  =  100,  the  effect  is  dramatic.  Both  extracted  velocity  com¬ 
ponents  are  significantly  less  accurate,  and  it  is  expected  that 


FIG.  12.  x=1.8,  Re=2.  <t=7,  b=2,  2=160.  Al=48.  Aj  =  10'^  /,=215,  2, =210,  Af  =  l.  (a)(i)  x  component;  (a)(ii)  y  component  for  J7=10'\  (b)(i)  x 
component;  (b)(ii)  y  component  for  17=10"‘;  (c)(i)  x  component;  (c)(ii)  y  component  for  »7=10"*. 
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FIG.  13  Here  ;r=1.8.  Re=2,  cr=7,  b  =  2,  7  =  160,  A/=48,  A^  =  10■^ 
;,  =  400.  7,  =  68,  and  M=  1.  (a)  Jf  component  of  velocity  for  7?=  10'*;  (b)  x 
component  of  velocity  for  77=10 


FIG.  15.  Here  A-=  1.8.  Re=2,  <r=7.  6  =  2.  7=160,  ^=48,  A^=10■^ 
7^=400,  and  7}-  10“*’.  Averaging  takes  place  on  the  M  temperature 

fields,  then  u  and  i;  are  extracted,  (a)  M=l;  (b)  M  =  100. 


for  larger  M,  the  results  will  be  further  degraded.  This  is 
because  noise  that  corrupts  quantities  appearing  nonlinearly 
in  the  hyperbolic  equations  can  lead  to  large  errors  in  the 
velocity  extracted  from  individual  corrupted  scalar  fields. 
Thus,  the  effects  on  the  means  of  the  extracted  velocity  com¬ 
ponents  can  be  profound. 


FIG.  14.  Rc=2,  <7=7,  ^  =  2,  J=160,  A/=48,  A5  =  10“^  /,=400,  J,=68, 
and  7?=  10"*’.  Averaging  takes  place  on  the  M  velocity  fields  extracted  from 
each  of  M  noisy  temperature  fields,  (a)  M  =  l;  (b)  M  =  100. 


A  better  method  for  ameliorating  the  effects  of  noise  is 
to  average  M  noisy  temperature  fields  and  to  extract  u  and  v 
from  the  mean  temperature  field.  By  fixing  7^-400,  J^  =  68, 
A5  =  10"^  and  ?7=10"^,  and  setting  M  =  1  and  100,  we  ob¬ 
tain  the  results  shown  in  Figs.  15(a)  and  15(b),  respectively. 
For  M  =  1 ,  the  errors  in  both  velocity  components  are  similar 
to  the  noise-free  case.  When  M  is  10,  100,  1000,  or  10  000, 
the  results  are  indistinguishable  from  the  noise-free  case. 
Therefore,  the  noise  in  the  temperature  field  has  no  effect  on 
the  extracted  velocity  for  M^IO.  Fixing  7^=400,  7^  =  68, 
A5  =  10~^,  increasing  77  to  10“''*,  and  setting  M  to  10,  100, 
1000,  and  10  000,  we  obtain  the  results  shown  in  Figs,  16(a), 
16(b),  16(c),  and  16(d),  respectively.  For  A7  =  l  [Fig.  13(a)], 
the  errors  in  both  velocity  components  are  large.  As  M  is 
increased,  the  errors  decrease  until  for  M  =  1 0  000,  the  error 
due  to  noise  is  insignificant  compared  to  those  associated 
with  discretization  onto  a  grid.  Therefore,  the  error  due  to 
noise  can  be  reduced  by  averaging  a  sufficient  number  of 
temperature  fields  before  extracting  the  velocity. 

VI.  DISCUSSION  AND  CONCLUSIONS 

We  have  shown  that  it  is  possible  to  extract  both  velocity 
components  from  advected  scalar  data  for  a  two-dimensional 
steady,  solenoidal  flow  with  appropriate  scalar  boundary 
conditions.  The  proof-of-concept  study  reported  here  shows 
that  we  can  deal  with  singularity  in  the  equations  from  which 
the  velocity  is  extracted,  as  indicated  in  Ref.  5. 

From  a  numerical  standpoint,  the  integration  along  the 
characteristic  curves  is  stable,  and  a  relatively  large  integra¬ 
tion  step  size  (100-1  (XX)  steps  per  characteristic)  is  suffi¬ 
cient.  For  the  flow  and  boundary  conditions  considered, 
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FIG.  16.  Here  ;r=  1.8,  Re  =  2.  cr=7,  =  7-  160,  N=48.  ^s=  10'%  /,=400,  ^=68,  and  Averaging  lakes  place  on  the  M  temperature  fields,  then 

u  and  i;  are  extracted,  (a)  M  =  10;  (b)  M  — 100;  (c)  M  =  1000;  (d)  Af  ==  10  000. 


knowledge  of  the  scalar  on  a  400X68  grid  is  sufficient. 
Available  CCD  arrays  can  be  used  to  obtain  scalar  data  with 
this  resolution. 

We  have  introduced  multiplicative  noise  and  have 
shown  how  to  deal  with  it  by  averaging  noisy  scalar  fields. 
More  sophisticated  filtering  and  smoothing  techniques  may 
be  more  attractive  at  higher  noise  levels.  The  trade-off  be¬ 
tween  the  number  of  images  that  need  to  be  processed  (with 
simple  averaging),  the  computational  complexity  (of  concern 
as  the  sophistication  of  the  image  processing  increases),  and 
the  cost  of  special  purpose  hardware  (to  reduce  the  CPU  time 
associated  with  image  processing)  will  determine  the  best 
approach  in  a  given  situation.  We  conjecture  that  additive 
noise  will  give  qualitatively  similar  results  at  low  noise  lev¬ 
els. 

The  ultimate  goal  of  this  work  is  to  develop  a  real-time 
version  of  this  velocity  extraction  technique  for  use  in  flow 
control  and  state  estimation.  In  unsteady  three-dimensional 
flows,  the  issue  of  spatial  resolution  will  be  critical  in  deter¬ 
mining  temporal  resolution  when  simultaneous  volumetric 
data  cannot  be  acquired,  in  which  case  two-dimensional  im¬ 
ages  must  be  “scanned”  in  the  third  direction.  The  spatial 
and  temporal  resolution  requirements  will  thus  be  important 
in  determining  framing  and  processing  rate  requirements  for 
use  of  this  technique  in  real-time  applications. 
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Abstract.  We  present  a  globalization  of  the  Krawczyk  algorithm  to  compute  all  real 
isolated  solutions  of  systems  of  N  real  polynomial  equations.  This  is  accomplished  by 
transforming  the  original  system  to  an  augmented  system  in  in  which  each  of  the  first  N 
variables  lies  in  the  interval  [-1,1],  and  the  (A^+l)-th  variable  lies  in  [0,°°).  The  domain  of  this 
latter  variable  can  be  divided  into  two  intervals  [0,1]  and  [1,°°),  the  second  of  which  is  mapped 
to  [0,1].  Thus,  the  entire  domain  can  be  examined  by  considering  the  finite  domain  [-1,1]^ 

X  [0,1]  for  each  of  two  systems.  One  of  the  augmented  systems  can  have  one  or  more  solutions 
not  corresponding  to  finite  solutions  of  the  original  system.  We  use  techniques  from  algebraic 
geometry  to  transform  the  (7'/+l)-th  dimensional  system  so  that  spurious  solutions  are  excluded, 
thus  restricting  the  solutions  to  those  corresponding  to  finite  solutions  of  the  original  system. 
The  algorithm  requires  bounds  for  multivariate  polynomials  over  a  finite  domain.  The  best  of 
three  bounding  methods  considered  uses  an  interval  extension  of  the  system,  which  is  stored  as  a 
rooted,  ordered  tree,  and  is  equivalent  to  an  //-dimensional  Horner  scheme  that  takes  advantage 
of  polynomial  sparsity. 

We  extend  the  approach  to  deal  with  systems  having  nonsimple  solutions  (i.e.,  with 
multiplicity  greater  than  unity),  at  which  the  Jacobian  vanishes.  This  is  accomplished  by 
constructing  a  new  system  that  reduces  by  one  the  multiplicity  of  nonsimple  solutions.  Except  in 
degenerate  cases,  this  approach  can  be  applied  sequentially  until  no  roots  have  multiplicity 
exceeding  one. 
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1.  Introduction.  We  present  here  a  new  computational  technique  for  finding  all  real 
solutions  in  of  systems  of  polynomial  equations 
(1)  =  0  \<i<N  . 

To  put  the  present  work  in  context,  we  begin  with  some  background. 

Nonlinear  equation  systems  such  as  (1)  arise  in  many  applications.  For  example,  a 
fundamental  problem  in  computer-aided  geometric  design  is  the  efficient  computation  of  all  the 
real  solutions  of  a  system  of  multivariate  polynomial  equations  within  a  given  domain.  Such 
equation  systems  can  arise  from  the  need  to  find  all  points  of  intersection  between  pairs  of 
curves,  and  curves  and  surfaces,  as  well  as  points  lying  on  curves  corresponding  to  the 
intersection  of  two  surfaces  [1].  Nonlinear  equation  systems  also  arise  in  computer-aided 
process  design  where  large  equation  systems  are  used  in,  for  example,  steady-state  process 
flowsheet  modeling  [2-3].  Another  important  application  is  the  solution  of  nonlinear  systems 
arising  from  spatial  discretization  of  boundary  value  problems  having  polynomial  nonlinearities, 
such  as  the  Navier-Stokes  equations. 

Solution  techniques  for  nonlinear  equation  systems  divide  into  two  broad  classes.  The  first 
is  the  class  of  local  methods,  where  one  seeks  a  solution  by  starting  with  an  initial  iterate.  These 
methods  include  single-  and  multi-point  iteration  techniques  (e.g.,  Newton-Raphson,  Muller's 
method).  The  second  class  of  methods  systematically  attempts  to  find  all  solutions  in  a  domain. 
These  methods  include  homotopy  and  subdivision  (e.g.,  bisection)  techniques,  as  well  as 
techniques  that  reduce  the  dimension  N  to  unity  (e.g.,  methods  using  resultants,  Grdbner  bases). 
If  enough  information  is  available  to  select  a  good  initial  iterate,  single  point  iteration  schemes 
can  find  single  solutions  well.  If  the  initial  iterate  is  vydthin  the  region  of  convergence  of  (1),  the 
method  will  converge.  If,  on  the  other  hand  the  initial  iterate  is  outside  the  region  of 
convergence,  the  method  might  converge  to  another  solution  or  not  converge  at  all. 

Algebraic  schemes  can  be  used  to  reduce  solution  of  (1)  to  finding  zeros  of  a  univariate 
polynomial.  The  main  approaches  use  elimination  [4-5]  or  the  Grobner  transTormation  [6-7]. 
The  advantage  of  these  methods  is  that  computation  of  the  solutions  of  a  univariate  polynomial 
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equation  is  a  standard  procedure,  so  that  these  techniques  are  capable,  in  principle,  of  finding  all 
solutions  in  .  The  manipulative  burden  of  reducing  an  //-dimensional  system  to  a  lower¬ 
dimensional  system  (typically  one  dimension)  can  be  significantly  diminished  by  symbolic 
mathematical  programs.  The  disadvantages  are  that  these  methods  suffer  from  numerical 
instability  arid  undesirable  cost  grovrth  (e.g.,  memory  and  time)  as  N  increases. 

Homotopy  methods  [8]  embed  (1)  in  a  family  of  systems  A  frequently  used 

embedding  is  the  simple  one-parameter  linear  homotopy  Kix,t)  =  (\-t)G{x)  +  tf{x).  The 
continuation  parameter  t  connects  the  target  system  (1)  to  a  system  =  for  which  all 
solutions  in  are  known.  An  advantage  of  homotopy  approaches  is  that  they  allow  all 
solutions  (real  and  complex)  of  (1)  to  be  found.  However,  since  the  starting  points  of  real 
solutions  cannot  be  identified  in  advance,  all  solution  paths  must  be  followed  in  order  to  compute 
all  the  real  solutions.  For  polynomial  systems,  the  number  of  paths  is  typically  equal  to  the  total 
degree. 

The  class  of  global  subdivision  techniques  includes  interval  methods.  These  methods  find 
all  solutions  of  (1)  in  a  finite  domain.  As  discussed  by  Moore  and  Jones  ([9],  p.  1056),  three 
possibilities  arise  when  examining  such  a  finite  domain.  One  is  sometimes  able  to  show  that  no 
solution  exists,  or  that  at  least  one  solution  exists.  Tbe  third  possibility  is  that  no  such  conclusion 
is  possible.  In  the  latter  case,  bisection  is  used  to  subdivide  the  given  domain  until  each 
subdomain  is  shown  to  either  contain  or  not  contain  a  solution.  A  key  advantage  of  this 
technique  is  that  it  finds  all  real  solutions  in  a  finite  domain  without  computing  the  complex 
solutions.  Interval  methods  that  take  advantage  of  large-grained  parallelism  now  provide  an 
efficient  approach  to  compute  all  solutions  in  a  finite  domain  for  application  problems  of  small 
and  intermediate  dimension  (cf.  [10-13]). 

In  this  paper,  we  show  how  a  technique  that  uses  interval  arithmetic  and  a  contraction 
mapping  algorithm  [9,14-15]  to  find  all  solutions  of  (1)  in  a  finite  domain  can  be  generalized  to 
find  all  solutions  of  (1)  in  .  The  generalization  to  unbounded  domains  is  significant  in  that  it 
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lifts  the  requirement  that  bounds  be  established  for  the  solution(s),  frequently  by  techniques  that 
are  ad  hoc  or  problem-dependent. 

2.  Development  of  a  global  method.  In  this  section,  we  introduce  notation,  review  a  few 
fundamentals  of  interval  analysis,  and  introduce  the  Krawczyk  operator.  We  then  proceed  to 
globalize  the  latter,  present  the  basic  algorithm,  and  discuss  how  it  deals  with  spurious  solutions, 
solutions  on  or  near  the  boundaries  of  intervals,  and  nonsimple  solutions. 

2.1.  Background.  Here,  we  establish  notation  and  briefly  review  the  basics  of  interval 
analysis  and  computation  [16-17],  the  Krawczyk  algorithm  [17-18],  and  Moore's  approach  to 
computing  real  solutions  of  multivariate  equations  by  bisection  of  interval  vectors  [9,17]. 

An  interval  A'  is  a  closed  and  bounded  set  of  real  numbers  [a,b].  An  interval  vector  2L  is 
an  ordered  «-tuple  of  intervals.  (In  this  paper,  scalars  and  vectors  are  denoted  by  Roman 
variables  with  no  underscoring  and  single  underscoring,  respectively.  Real  matrices  are 
represented  by  double-underscored  Greek  variables,  with  the  exception  of  the  N  x  N  identity 
matrix,  denoted  by  7^.  All  other  double-underscored  Roman  variables  denote  interval 
matrices.)  Arithmetic  operations  on  intervals  are  defined  straightforwardly.  For  example,  the 
product  Z  =  XY  defines  an  interval  containing  the  product  xy  for  all  x  e  A'  and  all  y  6  F .  An 
interval  extension  7^  of/is  an  interval-valued  function  whose  range  F{X)  includes  /(x)  for  all 
X  G  One  way  to  obtain  F  is  to  replace  the  real  variables  in /by  their  corresponding  interval 
variables.  Such  an  F  might  also  contain  points  that  do  not  correspond  to  /(x)  for  any  x€2L- 
We  denote  the  width  of  an  interval  X  by  w{X)  =  b  -  a.  The  midpoint  of  an  interval  is  defined 
as  miX)  =  {b  +  a)l2. 

The  algorithm  used  by  Moore  and  Jones  [9]  is  based  on  the  work  of  Krawczyk  [18],  who 
used  the  Brouwer  fixed  point  theorem  to  develop  an  interval-based  Newton-Kantorovich 
mapping.  The  Krawczyk  operator  is  defined  by 

(2)  K{JO  =  h-0fih)  +  {l^-0JiX)}{X-h)  , 

where  0  is  an  arbitrary  nonsingular  real  matrix  and  h  is  the  real  vector  m(20,  while  7^  is  the 
identity  matrix  and  £{X_)  is  an  interval  extension  of  the  Jacobian  of  the  vector-valued  function 
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/.  (Midpoints  of  interval  vectors  and  matrices  are  defined  as  the  midpoints  of  their  elements.) 

We  define  v(2C)  =  11^/  volume  of  a  rectangular  portion  of  . 

1=1 

2.2.  Globalization.  In  order  to  globalize  the  Moore-Krawczyk  algorithm  [9],  one  needs  to 
be  able  to  examine  the  infinite  domain  .  With  this  in  mind  the  polynomial  system  (1)  is 
mapped  from  the  original  system  (1)  in  to  an  augmented  system  in  R^*'  by  the 
transformation 


(3a) 


1-1/2 


L/=1  J 

(3b)  =  . 

After  multiplying  each  transformed  equation  by  where  d,  is  the  degree  of  the  /-th  equation, 
we  get 


(4a)  =  ^  , 

in  which  the  term  or  terms  of  highest  degree  in  each  equation  of  system  (1)  are  invariant.  The 
(A+l)-th  equation  is 

(4b)  g//+i(yp"'5y/,F)  =  XTi  ~  ^  • 

1=1 

Thus,  the  solutions  of  (4a, b)  lie  on  a  semi-infinite  right  circular  cylinder  of  unit  radius  for  N  =  2, 
and  on  a  semi-infinite  unit  hypercylinder  for  N>3.  The  rectangular  domain  to  be  examined  in 
order  to  find  all  solutions  of  (4a,b)  is 

(5)  T,  e[-l,l]  ,  e[0,oo)  ,  \<i<N  . 

Thus,  the  transformation  reduces  the  problem  to  one  for  which  only  one  variable  is  unbounded. 
We  note  that  if  a  null  solution  of  (1)  exists,  it  will  not  be  a  solution  of  (4).  The  existence  of  such 
a  solution  is  easily  verified. 

The  unbounded  domain  [0,«>)  can  be  broken  up  into  the  intervals  [0,1]  and  [1,°°).  Using 
the  transformation 


(6) 


w,  =  >/;  \<i<N 

^//+1  ~  V TW+I  9 
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the  interval  (1;°°)  is  mapped  to  [0,1],  This  mapping  gives  rise  to  the  system 

(7)  =  0  \<j<N+]  . 

Thus  the  entire  domain  can  be  examined  by  considering  a  finite  domain  for  each  of 
the  systems  (4a, b)  and  (7). 

2.3.  Basic  Algorithm.  An  algorithm  implementing  these  ideas  is  represented  by 
pseudocode  in  Figure  1.  The  algorithm  examines  finite  domains  7,  and  eliminates  those  that  can 
be  shown  to  not  contain  a  solution.  If  the  mapping  (2)  gives  a  contraction  {KiY)  C  7),  then  7 
contains  a  solution.  In  the  event  that  no  existence  or  nonexistence  result  is  obtained,  7  is 
bisected,  as  described  in  detail  below,  and  added  to  a  "stack"  of  finite  domains  to  be  examined. 
This  process  proceeds  until  the  initial  domain  has  been  examined  and  all  solutions  have  been 
found. 

Examination  of  each  finite  domain  begins  with  a  nonexistence  test  [9],  checking  to  see  if 

(8)  OgG/7) 

is  true  for  some  ;,  1  <  ;  <  A  + 1 .  If  (8)  is  satisfied  for  some  j,  then  7  contains  no  solutions  of 
(4a, b),  and  7  is  discarded.  If  (8)  is  false  for  each  j,  the  possibility  remains  that  7  contains  a 
solution  of  (4a,b).  In  that  case,  a  second,  new  nonexistence  test  (see  Appendix)  checks 

(9)  1«,(»>(D)1>  X  O  + 

for  each  j.  If  (9)  is  false  for  each  j,  we  compute  the  mapping  (2),  using  the  choice 

[/n(7(7*))]'’  if  f/,,  -^*7(7!)|<||/^  -^*^-^7(7*-’)|| 
otherwise 

Otherwise,  7  contains  no  solutions  and  is  discarded.  If  m{£)  is  singular,  we  bisect  the 
corresponding  7,  and  place  the  resulting  interval  vectors  on  the  stack.  Otherwise,  (2)  is 
computed  and  a  third  nonexistence  test  is  applied,  which  checks  for  emptiness  of 
(11)  Z  =  KiY)nY  . 

'  If  (1 1)  is  empty,  then  7  contains  no  solutions  of  (4a,b)  and  is  discarded.  If  (1 1)  4s  nonempty,  we 


check  to  see  if 


(12)  YjQKjH) 

is  true  for  1  <  ;  <  7/  + 1 .  If  so,  7  is  bisected.  Otherwise,  an  existence  test  which  checks 

(13)  Kj(Y)czY^ 

is  applied.  If  (13)  is  true  for  1  <  j<N  +],  a  contraction  has  taken  place  and  Y  is  guaranteed  to 

contain  a  solution  of  (4a, b).  In  that  event,  any  point  in  Y  will  be  a  safe  starting  point  (with 

guaranteed  convergence)  for  use  with  interval  Newton-Raphson  [9].  Details  of  the  convergence 

and  solution  acceptance  are  shown  in  Figure  1. 

At  this  point,  either  a)  a  solution  has  been  shown  to  exist  in  the  current  interval  vector  7,  b) 

nonexistence  tests  have  excluded  a  solution  in  7,  or  c)  nonexistence  tests  have  failed  to  exclude 

all  or  part  of  7  from  further  consideration.  The  final  step  is  deciding  whether  to  place  (11)  on 

the  stack,  or  bisect  (11)  and  place  the  resulting  interval  vectors  on  the  stack.  We  decide  this  by 

computing  v(Z)/v(7).  If  this  volume  ratio  is  less  than  y  (chosen  arbitrarily),  (1 1)  will  be  placed 

on  the  stack,  and  bisected  otherwise.  This  ends  a  typical  pass  through  the  main  algorithm. 

(Trial-and-error  experimentation  shows  that  the  threshold  y  - 1/2  works  well.) 

In  order  to  bisect  7,  the  direction  of  bisection  must  be  chosen.  This  choice  can  greatly 

affect  the  efficiency  of  the  overall  scheme.  We  consider  two  schemes  for  making  this  decision. 

Scheme  1  bisects  in  the  direction  with  index ]  for  which  >Vy(7^)  is  maximum.  This  scheme  tends 

to  minimize  the  ratio  max  w,(7,)/ min  wAYj).  Scheme  2,  described  in  [19],  defines 

j  '  j 

(14)  5^  =  max||7;j*|  ,  |7„j;fc|}w(7y)  , 

and  bisects  in  the  direction  with  index  j  for  which  Sj  is  maximum,  where  7/^^  and  J^  jj^  are  the 
lower  and  upper  bounds  of  the  jk-th  element  of  the  interval  Jacobian,  respectively.  A 
comparison  of  results  for  these  two  bisection  schemes  will  be  presented  in  section  4.4. 

2.4.  "Spurious"  Solutions.  In  some  cases,  some  solutions  of  the  augmented  system  (4a,b) 
will  not  correspond  to  finite  solutions  of  the  original  system  (1).  These  spurious  solutions,  with 
7w+i  =  Oj  are  introduced  by  the  transformation  and  multiplication  described  in-section  2.2.  For 
example,  the  system 
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(15a) 

(15b) 

(15c) 


5xf  -  6x1x1  +  x-^xl  +  2x,X3  =  0 
-2x^X2  +  2x^x1  +  2^3X3  =  0 

xf  +  x\  -  0.265  =  0 


when  transformed  gives  the  augmented  system 


(16a) 

(16b) 

(16c) 

(16d) 


5>^’  -  6ylyly]  +  y^y'y'  +  2y^y2yl  =  0 

-2y^y2 + 2y^y^2yl + ^yiy^yl  =  0 

>;f  +  y3'-0.265>^,'  =  0 

y^  +  yl  +  yl-'^  =  ^  > 


which  has  two  spurious  solutions  (0,0, ±1,0),  corresponding  to  unbounded  solutions  of  the 
original  system  (1 5a-c).  On  the  other  hand,  the  system 


(17a)'  1  +  2x,  -  3xf  +  4x2  +  ^^2  ~  6x,X2  =  0 

(17b)  -7-8x,  +  9x1^  -lOxj  -llxj  +12x,X2  =0 


when  transformed  gives  the  augmented  system 
( 1 8a)  yl  +  2y^y2  -  +  4y2y2  +  Sy^  -  6y^y2  =  0 

( 1 8b)  -7>'3  -  8y,>^3  +  9yf  - 1  Oy^j^^a  “  1 1>'2  + 1  =  0 

(18c)  >;,'+>/2'-l  =  0  . 


The  transformed  system  (18a-c)  has  no  real  solutions  with  >^3  =  0,  so  globalization  introduces  no 
spurious  solutions. 

It  is  apparent  from  the  properties  of  the  mapping  that  every  real  finite  solution  of  (1)  is  also 
a  real  finite  solution  of  (4a,b).  Thus,  the  real  solutions  of  (4a, b)  consist  of  the  real  solutions  of 
(1)  and  any  spurious  solutions  introduced  by  the  transformation  and  multiplication  described  in 
section  2.2.  No  such  solutions  occur  if  the  equations  formed  by  omitting  all  but  the  leading- 
order  terms  (i.e.,  those  of  highest  total  degree  [7])  in  (1)  have  only  isolated  solutions,  which 

N 

when  transformed  do  not  satisfy  =1.  Spurious  solutions  will  occur  if  at  least  one  original 

1=1 

variable  (x,,  -,X/^)  is  absent  from  the  leading-order  terms  in  each  equation  in  (1),  or  appears  in 
all  of  the  leading-order  terms  of  more  than  one  equation.  In  the  example  (15a-c),  X3  is  absent 
from  the  leading-order  terms  in  each  equation,  whereas  each  variable  appears  in  at  least  one 
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leading-order  term  in  at  least  one  equation  of  (17a,b),  in  which  system  no  variable  is  present  in 

all  leading-order  terms  in  more  than  one  equation. 

It  is  possible  to  transform  (4a, b)  so  that  solutions  with  =  0  are  excluded,  thus 

restricting  the  solutions  to  those  corresponding  to  finite  solutions  of  (1).  This  is  accomplished 

using  techniques  from  algebraic  geometry.  To  remove  all  solutions  with  =  0  from  (4a, b), 

the  ideal  quotient  [7]  of  the  ideals  generated  by  yfj^^  =  0  and  the  N  equations  of  (4a)  is  found. 

The  ideal  quotient  is  itself  an  ideal  defining  a  new  system  of  N  equations.  That  system, 

combined  with  (4b),  has  a  solution  set  identical  to  the  solution  set  of  (4a,b),  less  the  spurious 

solution(s)  of  the  latter.  To  compute  the  ideal  quotient  one  determines  the  intersection  of  two 

ideals,  using  Grobner  bases  [7].  For  example,  the  system  (15a-c)  has  twelve  real  solutions  in  . 

The  solution  set  of  the  augmented  system  (16a-d)  consists  of  fourteen  real  solutions,  including 

the  two  spurious  solutions  noted  above.  These  spurious  solutions  can  be  excluded  using  the 

above  technique.  The  new  system  is 

(19a)  68719476736-10V1-20722660655104-]0V3>'4- 

588558700544-10^3:^4  -3982698885>^,^  =0 

(19b)  2048000>'3^  -  (>Al>Q12y^y,  -  9511  y\  =  0 

( 1 9c)  64>',^  -t-  64>’j  - 1 7>^4  =  0 

(19d)  >^?+:V:+>'3-l  =  0  • 

The  solution  set  of  (19a-d)  consists  of  all  nonspurious  solutions  of  (16a-d),  i.e.,  those  with 
>'4  0.  The  main  effort  in  computing  an  ideal  quotient  is  computing  a  Grobner  basis.  As  noted 

in  section  1,  the  Grobner  transformation  can  be  used  to  find  all  real  solutions  of  (1).  An 
important  point  to  note  is  that  to  use  the  Grobner  transformation  a  specific  term  ordering  (e.g., 
>  XjXj  using  total  degree  ordering)  must  be  chosen  for  the  monomials  [7]  in  each 
polynomial  equation  of  (4a,b).  The  choice  of  a  term  ordering  (cf.  [6])  affects  the  computational 
complexity  of  constructing  a  Grobner  basis.  Our  experience  also  shows  that  the  ordering  of  the 
variables  (corresponding  to  permutations  such  as  (XpXjjXj),  (XpXjjXj),  etc.)  can  also 
significantly  influence  the  operation  count.  In  most  cases,  the  lexicographical  term  ordering 
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needed  to  compute  a  reduced  Grobner  basis  (most  useful  in  computing  multivariate  polynomial 
zeros)  requires  many  more  operations  than  to  compute  a  generic  Grobner  basis  using  some  other 
term  ordering.  Thus,  in  our  removal  of  spurious  solutions  from  (4a, b),  we  will  use  a 
nonlexicographical  term  ordering,  and  avoid  computing  any  reduced  Grobner  basis. 

In  closing  this  section  we  note  that  ours  is  a  "black  box"  algorithm  for  general  use.  On  the 
other  hand,  one  can  use  our  approach  as  a  "tool",  without  removing  spurious  solutions  from 
(4a,b).  In  that  event,  there  can  occur  cases  in  which  a  large  number  of  very  small  interval  vectors 
(with  width  less  than  a  tolerance  set  by  the  user)  contain  spurious  solutions.  It  has  been  our 
experience  that  these  interval  vectors  can  be  safely  discarded  without  loss  of  finite  real  simple 
solutions  (i.e.,  those  with  unit  multiplicity). 

2.5.  Multiply-Computed  and  Boundary  Solutions.  This  contraction-mapping  algorithm 
will  not  find  a  solution  lying  on  the  boundary  of  the  initial  domain  or  on  a  boundary  created  by 
bisection.  To  ensure  that  all  solutions  of  (1)  are  found,  we  take  two  precautions  to  include  the 
boundaries  of  each  domain.  First,  when  the  interval  vector  Y  with  components  Y,  =  /,>'/,„]  is 

bisected  in  the y-th  direction  we  do  so  according  to 

(20a)  F,,  =  and  Y:^^  =  [y,„y,J  \<i<N+\  ,  i^j 

(20b)  Y,j  =  [yjj,  (1  /  2  -  £, )  •  y, ,  +  (1  /  2  +  e, )  •  J 

(20c)  Y^j  =  [(1  /  2  +  £,)  •  yjj  +  (1  /  2  -  £,)  •  y,. J  , 

so  that  no  point  originally  in  the  interior  of  F  lies  on  the  boundary  of  either  new  subdomain. 
Here,  £,  >  0  is  typically  0.005.  The  second  step  involves  expansion  of  the  interval  vectors  F, 
and  Kj  before  they  are  placed  on  the  stack.  The  new  interval  vectors  are 
(2 1  a)  F*,  =  [y,jj  -  £2^(7, ,),  y, ,  „  +  £2w(F, _,)] 

1  <  I  <  JV  + 1 


(2 1  b)  F2* ,  =  [^2.,.,  -  y2.i.u  +  £2>^(F2./)] 

where  £2  >  0  is  typically  0.0025.  Thus,  any  solution  lying  on  a  boundary  of  F,  or  F2 

•  • - ♦  — 

included  in  F,  or  £2 . 
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These  "precautions"  can,  however,  lead  to  computing  a  simple  solution  more  than  once. 
This  problem  is  addressed  after  all  solutions  of  (1)  have  been  found.  We  first  compute  the 
maximum  distance  betw'een  all  pairs  of  solutions.  Then  if  the  distance  dj  between  a  given 
pair  of  solutions  satisfies 

(22)  ^  ^3  ’ 

we  consider  the  possibility  that  they  are  the  same  solution,  where  a  typical  value  of  £3  is  10“*. 
This  list  can  contain  simple  solutions  computed  more  than  once.  We  place  a  box  around  each 
cluster  of  such  solutions  and  compute  the  solution(s)  in  this  box  using  the  algorithm  discussed  in 
sections  2. 1-2. 3.  Solutions  which  are  unique  will  be  thus  identified.  The  final  solution  set 
contains  the  original  unclustered  solutions  and  those  solutions  found  by  this  "declustering" 
process. 

2.6.  Nonsimple  Solutions.  Like  every  Jacobi an-based  procedure,  the  method  described 
above  encounters  difficulties  for  nonsimple  solutions  (i.e.,  those  with  multiplicity  greater  than 
one),  at  which  points  the  Jacobian  is  singular.  Since  the  determinant  of  the  Jacobian  is  a  single 
function  of  the  variables  and  vanishes  at  every  nonsimple  solution,  it  is  easy  to  see  that  at  each 
such  point  the  Jacobian  either  has  an  isolated  singularity  or  is  singular  on  one  or  more  curves  or 
surfaces  passing  through  the  point.  Thus,  unless  the  Jacobian  has  only  isolated  singularities, 
Newlon  iteration  will  encounter  curves  (in  the  bivariate  case)  or  surfaces  or  hypersurfaces  (more 
generally)  on  which  the  Jacobian  is  singular.  This  typically  leads  to  poor  convergence  behavior 
[20-21], 

Here  we  show  how  to  construct  a  new  system  whose  solutions  include  all  of  the  nonsimple 
solutions  of  the  original  system.  The  multiplicity  of  these  solutions  is  typically  reduced  by  one. 
If  no  solution  of  the  original  system  has  multiplicity  greater  than  two,  then  the  problem  is 
reduced  to  one  that  can  be  dealt  with  using  the  techniques  described  in  sections  2. 1-2. 5. 
Otherwise,  the  reduction  described  here  can  be  applied  sequentially  until  the  multiplicity  of  the 
solution  with  highest  original  multiplicity  has  been  reduced  to  one,  allowing  reliable  numerical 
solution.  Our  approach  is  a  direct  extension  of  a  technique  for  the  univariate  case  [22]. 
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We  construct  a  new  system  of  equations  consisting  of  some  ( N  -  l)-element  subset  of  the 
original  equations,  augmented  by  the  determinant  of  the  original  Jacobian 
(23)  /,„=U(/„/,.-./»)l=0  . 

As  indicated  in  examples  below  (section  4.4),  this  process  can  be  repeated  to  handle  nonsimple 
solutions  with  multiplicity  greater  than  two.  We  note  here  that  the  new  system  can  have  one  or 
more  solutions  not  satisiying  (1).  Also,  there  are  N  ways  to  construct  the  new  equation  set.  As 
discussed  in  example  5  below,  some  are  more  useful  than  others. 

The  method  performs  well  for  almost  all  systems  with  nonsimple  solutions,  as  indicated  by 
the  examples  below.  As  discussed  at  the  end  of  this  subsection,  it  does  not  work  for  Powell's 
degenerate  system  (cf.  [23]),  which  arises  from  applying  derivative-based  techniques  to 
minimization  of  "Powell's  singular  function"  ([24],  p.  150),  and  another  example  we  have 
developed. 

1.  A  bivariate  system  with  a  solution  of  multiplicity  two: 

(24a)  /,  =x,^ +  X2 -2  =  0 

(24b)  /z  =  +  lx]  -h  3x,  -H  Xj  -  7  =  0  . 

This  system  has  one  distinct  real  solution  (1,1),  of  multiplicity  two. 

The  reduced-multiplicity  system  {f2,ff)  is 
(25a)  x]  +  2x2  +  3x,  -t-  Xj  -  7  =  0 

(25b)  4x,X2-6x2+2x,  =0  , 

with  a  real  simple  solution  (1,1).  It  is  easily  shown  that  the  other  real  solution  of  (25a, b)  does 
not  satisfy  (24a,b).  In  this  example  and  the  next  three,  any  choice  of  - 1  equations,  augmented 
by  the  Jacobian  of  the  ^/-dimensional  system,  reduces  the  multiplicity. 

2.  A  bivariate  system  with  a  solution  of  multiplicity  three: 

(26a)  /,  =  xj*  -I-  X2  -  2  =  0 

(26b)  /2  =  +  2x1  +  3  ^  • 

This  system  has  two  distinct  real  solutions  (1,1)  arid  (-1,-1),  the  latter  having  muRiplicity  three. 
We  construct  the  new  system  in  two  steps.  In  the  first,  we  compute 
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(27a)  /3=|^(/iJ2)|  =  4^i^2- 2^1 -2^2=0  , 

and  in  the  second  step,  for  the  pair  (/jj/s),  compute 

(27b)  /4=|^(/iJ3)l  =  8^f-8^2-4^i+4A:2=0  . 

Then  the  new  system  (/i ,  /4)  is 

(28a)  x^  +  xl-2  =  0 

(28b)  8xi* -8x2^ -4j:, +4x2  =0  • 

The  new  system  (28a,b)  has  the  same  solutions  as  (26a, b),  with  each  being  simple. 

3.  A  trivariate  system  with  a  solution  of  multiplicity  two: 

(29a)  /,  =  x,^  +  X2  -  2  =  0 

(29b)  f2  =  +  2X2  +  3x,  +  Xj  -  7  =  0 

(29c)  /3  =  xi'-X2^  +  X3+2  =  0  . 

The  new  system  is 

(30a)  xj*  +  2X2  +  3x,  +  X2  -  7  =  0 

(30b)  xj^-Xj +X3+2  =  0 

(30c)  ■  2x,  -  6x2  +  4x,X2  =  0  . 

Here,  the  original  system  has  one  distinct  real  solution  of  multiplicity  two  at  (1,1, -2).  We  find 
this  as  a  simple  solution  of  the  new  system  (30a-c). 

4.  A  bivariate  system  with  two  solutions  of  multiplicity  two: 

(31a)  /,  =  xf +  X2 +x, -2x2 -5  =  0 

(31b)  /2  =  52xf +73x2'  + 72x,X2  -20x,- 110x2 -575  =  0  . 

The  new  system  is 

(32a)  xf  +  Xj  +  X,  -  2X2  -  5  =  0 

(32b)  -144x2^  +  144xf+84x,x2  +  60x,  + 3 30x2 -150  =  0  . 

Here,  the  original  system  has  two  distinct  real  solutions  at  (-2,-1)  and  (1,3),  each  of  multiplicity 
two.  We  find  these  as  simple  solutions  of  the  new  system  (32a,b). 

5.  A  bivariate  system  with  a  solution  of  multiplicity  two  f21J: 

(33a)  /j=x  +  /=o 
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(33b)  f2  =  i^+y^+y^=^  • 

The  new  system  (/ij/s)  is 

(34a)  x  +  y^=0 

(34b)  |x  +  2y  =  0  . 

The  original  system  has  two  distinct  real  solutions  at  (-4,2)  and  (0,0),  with  the  latter  having 
multiplicity  two.  The  simple  solution  at  (-4,2)  is  found  directly  from  (33a,b),  while  the  solution 
of  multiplicity  two  at  (0,0)  is  found  from  (34a,b).  In  this  example,  the  choice  {/],/■),)  leads  to 
the  desired  multiplicity  reduction,  while  {fj,  f 3)  does  not. 

6.  A  bivariate  system  with  a  solution  of  multiplicity  three  /20J: 

(35a)  /,  =  xf -i- =  0 

(35b)  f2=X2  +  xl=0  . 

This  system  has  four  distinct  real  solutions  at  (0,-1),  (-1,-1),  (-1,1),  and  (0,0),  with  the  latter 
having  multiplicity  three.  As  in  example  2,  one  stage  of  multiplicity  reduction  (using  /2  and 
f-^)  results  in  a  system  where  the  solution  of  (35a, b)  with  multiplicity  three  is  now  a  solution  of 
multiplicity  two.  The  second  stage  (using  /2  and  =|i(/2./3)l)  gives 
(36a)  X2  +  X2  =0 

(36b)  X)  +  4xiX2  +  4x]X2  =  0  . 

The  simple  solutions  of  (35a, b)  are  found  directly  from  that  system,  while  (36a,b)  provide  the 
solution  of  multiplicity  three  at  (0,0)  and  the  simple  solution  at  (0,-1).  In  this  example,  {f\,f3) 
gives  no  multiplicity  reduction  in  the  first  stage,  and  {f-^,  ff)  gives  none  in  the  second  stage. 

We  have  found  two  test  cases,  each  with  nonsimple  isolated  solutions,  that  do  not  yield  to 
our  approach.  One  is  Powell's  four-dimensional  system  [23-24],  for  which  the  Jacobian  vanishes 
on  two  hyperplanes,  on  each  of  which  one  of  the  original  equations  also  vanishes.  Kearfott's 
algorithm  [25]  also  experienced  difficulty  with  this  case.  A  second  example  is 
(37a)  x'l-xl=0 

■  (37b)  •  XiX2='0  , 
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for  which  the  Jacobian  vanishes  only  at  the  origin.  In  each  case,  application  of  our  approach 
without  prior  multiplicity  reduction  yields  small  regions  (containing  the  solution)  for  which  no 
definitive  statement  is  possible.  Use  of  multiplicity  reduction  for  Powell's  example  fails  to 
define  a  new  system  for  w'hich  solutions  are  isolated.  For  (37a, b),  the  solution  of  the  new  system 
does  not  have  lower  multiplicity. 

3.  Bounding  multivariate  polynomials.  Our  approach  requires  bounding  each  equation  of 
(4a,b)  over  a  finite  domain  Y.  Three  methods  are  discussed  and  implemented  for  bounding 
multivariate  polynomials  with  real  coefficients  over  a  finite  domain  7. 

In  the  first  method  [26],  one  transforms  each  polynomial  of  (4a,b)  or  (7)  into  generalized 
Bernstein  form,  computes  its  Bernstein  coefficients,  and  takes  the  minimum  and  maximum  as  the 
bounds  over  7.  This  method  gives  exact  bounds  if  and  only  if  the  original  polynomial  is 
monotonically  increasing  or  decreasing  on  the  domain  over  which  a  bound  is  required. 

The  second  method  [26-27]  is  based  on  the  mean  value  theorem  and  uses  function 
evaluations  and  a  Taylor-like  correction  term  to  bound  each  polynomial  of  (4a,b)  or  (7)  over  7. 
This  method  relies  in  part  on  polynomial  evaluations.  The  evaluations  are  computed  using  the 
tree  traversal  method  [28].  To  use  this  method  the  coefficients  of  (4a,b)  or  (7)  are  stored  as  a 
rooted,  ordered  tree  and  then  traversed  by  pre-order  traversal.  This  technique  is  equivalent  to  an 
A'-dimensional  Horner  scheme  and  takes  advantage  of  polynomial  sparsity. 

The  third  method  [17]  is  based  on  the  interval  extension  Gj{Y)  of  either  (4a, b)  or  (7).  An 
interval  extension  is  found  by  evaluating  the  polynomial  (real  variables  replaced  with  interval 
variables)  using  the  tree  traversal  technique  used  in  the  second  method. 

The  polynomial  bounding  method  was  chosen  based  on  extensive  tests  for  univariate  and 
bivariate  cases.  It  was  found  that  the  first  method  gives  the  tightest  bounds  for  a  given 
computational  cost.  The  second  is  comparable  to  the  first  in  time  but  easier  to  implement  in 
higher  dimensions.  The  third  method  is  much  faster,  but  provides  less  tight  bounds.  The  first 
method  is  the  hardest  to  implement  using  the'  tree  traversal  representation  "of  polynomials 
employed  here.  The  second  and  third  methods  are  easily  implemented  using  some  form  of  tree 
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traversal,  while  also  being  easily  extended  to  higher  dimensions.  After  running  examples 
utilizing  each  method  in  the  "solver,"  we  conclude  that  the  third  method  gives  the  best  overall 
performance  (i.e.,  lowest  total  time). 

These  techniques  give  bounds  that  include  the  true  range  of  gj.  Computational  burden  is 
associated  with  the  fact  that  these  bounds  are  not  "light."  The  extent  to  which  this  is  acceptable 
depends  on  the  relationship  between  bounding  costs  and  other  computational  costs.  For  this 
reason,  it  is  desirable  to  find  a  means  of  increasing  the  tightness  of  any  bounding  scheme.  This 
is  accomplished  by  bisecting  the  initial  domain  and  computing  bounds  for  each  resulting 
subdomain.  The  final  bounds  over  the  initial  domain  are  taken  to  be  the  minimum  and  maximum 
of  the  bounds  found  for  the  subdomains.  This  is  continued  until  the  bounds  attain  the  desired 
tightness  over  the  initial  domain. 

4.  Computational  results.  In  this  section,  we  present  results  for  four  groups  of  test 
problems,  along  with  a  comparison  of  bisection  techniques.  The  first  group  of  test  problems 
consists  of  what  have  become  standard  tests  for  nonlinear  equation  solvers.  The  second  and  third 
sets  consist  of  A^-dimensional  quadratic  systems  with  random  coefficients,  and  bivariate  M-th 
degree  systems  with  random  coefficients. 

4.1.  Standard  Test  Cases.  The  first  test  suite  is  a  subset  of  the  test  suite  from  [25],  and 
includes  problems  with  a  physical  basis  (e.g.,  a  combustion  problem)  and  made-up  test  cases. 
These  problems  have  been  used  as  test  cases  for  homotopy  methods,  as  well  as  for  simplicial  and 
interval-Newton  methods.  The  "Remarks"  are  generally  taken  or  paraphrased  from  [25]. 

7.  Brown's  almost  linear  system  (N=5): 

(38a)  /,  =  x,  +  ix^.-(iV-H)  =  0  \<i<N-\ 

7=1 

(38b)  /^=nx,-i  =  o 

7=1 

Remarks:  The  Jacobian  matrix  is  said  to  be  ill-conditioned  at  the  two  solutions  of 
smallest  magnitude  [25].  Thisiiinction  [29]  is  also  found  in  [23]- 
It  is  shown  in  [23]  that  (38a, b)  has  two  solutions  for  even  N  and  three  for  odd  N.  We  find  three. 
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2.  Combustion  chemistry  problem: 

=-1691  10^2X4  +  2.177 -10’ X2  +  0.55x,x,  +  0.45x,  - =  0 
/2  =  1 . 585  •  X2X4  +  4. 1 26  - 1 0“^  X1X3  -  8. 285  •  10^  X1X4  +  2. 284  •  1 X3X4 

- 1.91 8  •  10''x:3  +  48.4x4  -  27.73  =  0 

f,  =  xl-x^=0 

Remarks:  These  equations  arise  in  a  chemical  equilibrium  problem  [30],  The  problem 
has  been  solved  via  continuation  methods.  It  has  a  unique  solution  within 
the  nonnegative  unit  box,  but  has  other,  nonphysical  solutions  in  larger 
domains. 

We  find  two  real  solutions,  one  physical  and  one  nonphysical.  (Note  that  the  coefficient  of  the 
X3  term  in  /2  is  given  with  different  signs  in  [25]  and  [30].) 

3.  A  high-degree  polynomial  system: 


(40a) 

5x’  -  6xf  Xj  +  XjXj  +  2X|X3  =  0 

(40b) 

-2xfx2  +2x,^X2  +2X3X3  =  0 

(40c) 

x?+x^ -0.265625  =  0 

Remarks:  This  problem  has  12  real  solutions,  which  are  all  within  the  box  [-0.6,0. 6]  x 
[-0.6,0. 6]  X  [-5,5],  and  eight  other  finite  solutions.  The  total  degree  is  126. 
Thus,  this  system  causes  trouble  for  most  homotopy  continuation  methods. 

We  find  all  12  real  solutions. 

4.  Rosenbrock's function: 

(41a)  /,  =l-x,  =0 

(41b)  /2  =  10(x2-xi')  =  0 

Remarks:  This  problem  is  a  standard  test  case  for  minimization  algorithms  [31]  and  has 
also  been  considered  from  the  standpoint  of  a  nonlinear  equation  system  [32]. 
We  find  the  sole  solution  (1,1). 

5.  A  variable-dimension  system  of  quadratics: 


(39a) 

(39b) 

(39c) 

(39d) 
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(42a) 

(42b) 


\<i<N-] 


/,  =  (X,  -  0.1)^ +  x,,,- 0.1  =  0 

fs  =(X;,  -0.1)^+  X,- 0.1=0 

Remarks:  This  simple  system  of  quadratics  can  be  used  to  test  the  effects  of  increasing 
dimension  on  bisection  methods.  Since  its  Jacobian  matrix  is  sparse,  it  can 
also  be  used  to  debug  techniques  for  handling  structured  problems. 

We  find  the  two  real  solutions  for  W  =4. 

4.2.  N-Dimcnsional  Systems  of  Quadratics.  This  test  suite  consists  of  systems  of 
quadratics 

(43)  fi  =  a,  +  ^b,jXj  +  ^  \<i<N 

y=l  y=l  *=1 

with  varying  dimension  N.  The  coefficients  are  chosen  randomly  between  -10  and  10,  so  that 
unlike  example  5  of  section  4.1,  the  Jacobian  is  typically  not  sparse.  This  test  suite  can  also  be 
used  to  test  the  effects  (e.g.,  on  cost  growth)  of  increasing  dimension.  Table  1  shows  running¬ 
time  statistics  on  a  166MHz  Pentium  machine  for  10  test  systems  for  each  N.  Running  time 
depends  strongly  on  dimension.  The  large  standard  deviations  indicate  a  considerable  sensitivity 
with  respect  to  the  coefficients. 

4.3.  Pairs  of  M-th  Degree  Equations.  This  test  suite  consists  of  bivariate  systems  of 
various  degrees.  The  coefficients  are  again  chosen  randomly  to  be  between  -10  and  10.  Table  2 
show's  running-time  statistics  (same  machine  as  used  in  section  4.2)  with  an  ensemble  of  100 
bivariate  systems.  One  can  see  by  comparing  tables  1  and  2  that  cost  grows  less  rapidly  with 
degree  than  with  dimension,  as  expected.  We  also  note  that  the  ratio  of  mean  times 
<T  >f^l<T  >M-^\  increases  slightly  with  the  degree  M,  and  that  although  the  ratio  of  the 
standard  deviation  to  the  mean  time  decreases  slowly  (and  not  quite  monotonically)  with  degree, 
the  ratio  of  maximum  to  minimum  run  times  decreases  consistently  and  significantly  with 
degree. 

4.4.  Bisection  Scheme  Comparison.  In  this  section  we  present  timing  results  for  the  two 
bisection  schemes  discussed  in  section  2.3.  Table  3  shows  running-time  statistics  (same  machine 
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degree.  Scheme  2  is  slightly  faster.  For  small  degree,  some  variation  in  running  time  is  expected 
for  test  suites  having  identical  coefficients. 

5.  Discussion.  For  broad  classes  of  multivariate  polynomial  equation  systems,  the 
numerical  solution  technique  discussed  in  this  paper  has  several  advantages  over  other  global 
methods.  First,  our  technique  seeks  and  finds  only  real  solutions  of  (1),  unlike  homotopy 
methods  which  must  track  all  solutions  (i.e.,  complex  and  real  solutions).  Second,  our  method 
can  be  used  to  find  all  solutions  in  a  finite  subdomain  of  (e.g.,  when  the  solutions  must  all  be 
nonnegative,  as  when  computing  concentrations  in  chemically  reacting  systems,  allocating 
resources,  etc.)  without  dealing  with  the  rest  of  ,  thus  making  it  possible  to  omit  the 
globalization  transformation  section  2.2.  Third,  the  technique  can  be  tailored  to  quadratic 
systems  (for  which  the  Jacobian  is  linear  and  can  be  bounded  exactly)  to  significantly  reduce  run 
time.  Finally,  by  assigning  subdomains  to  all  available  processors,  the  code  can  easily  be 
parallelized  [12]. 

Our  approach  to  computing  nonsimple  solutions  also  has  advantages  over  existing  methods. 
First,  unlike  approaches  [20-21]  that  carefully  select  the  initial  iterate  to  "coax"  Newton  methods 
to  converge  at  singular  points  (with  generally  sublinear  convergence),  we  consider  a  sequence  of 
new  systems,  each  having  a  lower  multiplicity  at  these  singular  points  until  the  multiplicity  is 
reduced  to  one.  This  allows  our  interval  extension  of  Newton's  method  to  converge 
quadratically.  We  also  note  that  the  multiplicity  reduction  technique  can  be  used  as  a  pre¬ 
transformation  for  Newton-Raphson  and  other  single-point  iteration  techniques. 
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APPENDIX 


A  computationally  verifiable  sufficient  condition  is  given  for  nonexistence  (in  a  rectangular 
domain)  of  a  solution  to  a  system  of  polynomial  equations.  Using  this  test,  which  is  based  on  the 
Krawczyk  operator,  the  number  of  bisections  required  by  the  Moore-Krawczyk  algorithm  can  be 
significantly  reduced. 

We  call  the  reader's  attention  to  the  explicit  dependence  of  (2)  on  the  real  matrix  0,  which 
can  be  chosen  arbitrarily,  subject  only  to  We  note  here  that  the  theorems  [9,17] 

concerning  application  of  (2)  are  true  for  arbitrary  nonsingular  0.  In  what  follows,  we  exploit 
this  flexibility  to  develop  a  nonexistence  criterion. 

We  begin  by  recognizing  that  every  solution  of  (1)  in  X_  must  lie  in  for  every 

h  in  X  and  every  nonsingular  0.  As  we  have  fixed  h,  we  will  concentrate  on  choosing  0. 
We  rewrite  (2)  in  component  form  (using  index  notation)  as 

(Al)  K,  =  m{X,)-0Jjih)^{5„  -  ^^.(/«(J,*)  +  [-l,l]X^y*))}(^*  -/”(^*)) 
where  J  is  in  midpoint-halfwidth  form  [17].  We  denote  and  X,  by  [^,  A^,]  and  {^,XX 
respectively,  and  note  that  K{X_,h,0){^X_  will  be  empty  (and  hence  (1)  will  have  no  solutions 
in  X)  if 

(A2)  or  A:,  <^ 

for  some  combination  of  i  and  (nonsingular)  0.  We  are  thus  led  to  rewrite  (Al)  as 
(A3)  K,  =  miX,)-0Jj{h)^{5^,  -  0ym{Jj,)- ^>(,XJ;*)[-l,l]}iw(;ir*)[-l,l] 
from  which  it  follows  that 

(A4a)  ^  =  m{X,)-0ijf -  ^yW2(-''y*)|  +  |^y|iw(J^*)}^w(2f*) 

and 

(A4b)  J.  =  miX^)-0,jfj{h)  +  ^d,,  -  0ymiJj,)\  +  \0^J\^wiJj,)}\w{X,) 

From  (A4a,b)  we  have 

(A5a)  ^  >  m{Xi)- +|^y|Jw(^y*)|  +  |^y|2>v(*7y*)}^ti'(2r*) 

and 
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(A5b)  +  , 

from  which  follows  an  important  result. 

Theorem.  If 

(A6)  |//4)|>(M-'a)  +  |«(-'a)|)!'^<'^») 

is  true  for  some  j,  then  (1)  has  no  solutions  in  2L- 

Proof. 

From  (A2)  and  (A5b),  it  follows  that  if 

is  true  for  some  i  and  nonsingular  0,  then  ^  •  Therefore,  from  (A7), 


will  be  sufficient  to  ensure  emptiness.  Assume  (A6)  holds  for  some  j.  There  are  two  cases. 
First,  consider  fj>0.  Then  choose 


(A9) 


o’iX.) 


>0 


for  one  m  (we  take  m  =  j  for  simplicity),  0„j  =  0  for  j,  and  =  <5„  for  j  Thus,  (A8) 


is  satisfied  and 
(A  10) 


For  /^  <  0 ,  we  need  only  choose 

> _  _ 


>0 


for  one  m,  with  the  other  elements  of  0  chosen  as  before.  Thus,  if  (A6)  is  satisfied  for  somey, 
we  can  construct  a  nonsingular  0  such  that  (A8)  is  satisfied,  and  K{2^,h,0)r\K  will  be  empty, 
proving  the  theorem. 

A  simple  interpretation  of  this  theorem  is  possible.  For  iV  =  1  (a  one-dimensional  system), 
(1)  will  have  no  zeros  on  an  interval  [a,b)  if 


(All) 


b-a 

- max 

2  *e[a.i] 


the  meaning  of  which  is  obvious.  In  N  dimensions,  satisfaction  of  (A6)  means  that  the  y-th 
component  of  the  function  value  f{h)  is  too  large  for  f  j  to  vanish  in 


23 


Table  1 


Running-time  statistics  for  quadratic  systems  as  a  function  of  dimension  N. 


dimension,  N 

maximum 
time  (sec) 

minimum 
time  (sec) 

mean 
time  (sec) 

standard 

deviation 

2 

0.187 

0.032 

0.113 

0.047 

3 

2.875 

0.562 

1.728 

0.731 

4 

133.8 

31.44 

52.54 

29.74 

5 

1267 

420.5 

949.4 

289.8 

Table  2 

Running-time  statistics  for  pairs  of  higher-order  equations  as  a  function  of  degree  M. 


degree,  M 

T 

max 

maximum 
time  (sec) 

T 

min 

minimum 
time  (sec) 

<T> 
mean 
time  (sec) 

standard 

deviation 

2 

0.281 

0.031 

0.079 

0.039 

3 

0.938 

0.125 

0.295 

0.135 

4 

1.781 

0.344 

0.823 

0.289 

5 

3.844 

0.875 

1.916 

0.715 

Table  3 

Mean  running-times  for  pairs  of  higher-order  equations  using  two  bisection  schemes. 


degree, M 

Scheme  1 
mean  time  (sec) 

Scheme  2 
mean  time  (sec) 

2 

0.087 

0.081 

3 

0.295 

0.275 

4 

0.82~3 

0.773 

5 

1.916 

1.709 

while  stack  not  empty 

get  interval  vector  Y  from  stack 
compute  Gj(Y) 

if  (8)  false 

no  solution  in  y 
else  if  (9)  true  for  some  j 
no  solution  in  y 
else  if  0  singular 

bisect  y  and  place  its  halves  on  stack 
else  if^(Dny  =  0 
no  solution  in  y 
else  if  yc^(y) 

bisect  y  and  place  its  halves  on  stack 
else  if  ^(y)  c  Y 
solution  lies  in 

ifl|£®  =  /^-£^(Z)||<l 

=  y 

for  /:=1  ’.^iterations 
do 

compute  0*  =  0(y* ) 
compute 

if  l|£*|l/li*"^ll>l 

end  if 

compute  y*'*’'  =y*  n^(y*) 

while  ||y*+’  -  >1||  >  £ 

end  for 
end  if 

else 

compute  Z  =  K(Z)  Z 
if  v(Z)/v(y)<y 
place  Z  on  stack 
else 

bisect  Z  and  place  its  halves  on  stack 

end  if 
end  if 
end  if 
end  while 

Figure  1 .  Pseudocode  of  the  basic  algorithm. 
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